ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-2.0/src/UseTheForce/DUFF.cpp

Comparing trunk/OOPSE-2.0/src/UseTheForce/DUFF.cpp (file contents):
Revision 1650 by gezelter, Tue Oct 26 22:24:52 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC

#	Line 1 | Line 1
1	<	#include <stdlib.h>
2	<	#include <stdio.h>
3	<	#include <string.h>
1	>	/*
2	>	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	>	*
4	>	* The University of Notre Dame grants you ("Licensee") a
5	>	* non-exclusive, royalty free, license to use, modify and
6	>	* redistribute this software in source and binary code form, provided
7	>	* that the following conditions are met:
8	>	*
9	>	* 1. Acknowledgement of the program authors must be made in any
10	>	*    publication of scientific results based in part on use of the
11	>	*    program.  An acceptable form of acknowledgement is citation of
12	>	*    the article in which the program was described (Matthew
13	>	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	>	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	>	*    Parallel Simulation Engine for Molecular Dynamics,"
16	>	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	>	*
18	>	* 2. Redistributions of source code must retain the above copyright
19	>	*    notice, this list of conditions and the following disclaimer.
20	>	*
21	>	* 3. Redistributions in binary form must reproduce the above copyright
22	>	*    notice, this list of conditions and the following disclaimer in the
23	>	*    documentation and/or other materials provided with the
24	>	*    distribution.
25	>	*
26	>	* This software is provided "AS IS," without a warranty of any
27	>	* kind. All express or implied conditions, representations and
28	>	* warranties, including any implied warranty of merchantability,
29	>	* fitness for a particular purpose or non-infringement, are hereby
30	>	* excluded.  The University of Notre Dame and its licensors shall not
31	>	* be liable for any damages suffered by licensee as a result of
32	>	* using, modifying or distributing the software or its
33	>	* derivatives. In no event will the University of Notre Dame or its
34	>	* licensors be liable for any lost revenue, profit or data, or for
35	>	* direct, indirect, special, consequential, incidental or punitive
36	>	* damages, however caused and regardless of the theory of liability,
37	>	* arising out of the use of or inability to use software, even if the
38	>	* University of Notre Dame has been advised of the possibility of
39	>	* such damages.
40	>	*/
41
5	–	#include <iostream>
6	–	using namespace std;
42
43	<	#include "UseTheForce/ForceFields.hpp"
9	<	#include "primitives/SRI.hpp"
10	<	#include "utils/simError.h"
11	<	#include "types/DirectionalAtomType.hpp"
43	>	#include "UseTheForce/DUFF.hpp"
44		#include "UseTheForce/DarkSide/lj_interface.h"
45	+	#include "UseTheForce/DarkSide/charge_interface.h"
46		#include "UseTheForce/DarkSide/dipole_interface.h"
47		#include "UseTheForce/DarkSide/sticky_interface.h"
48	+	#include "UseTheForce/ForceFieldFactory.hpp"
49	+	#include "io/DirectionalAtomTypesSectionParser.hpp"
50	+	#include "io/AtomTypesSectionParser.hpp"
51	+	#include "io/LennardJonesAtomTypesSectionParser.hpp"
52	+	#include "io/ElectrostaticAtomTypesSectionParser.hpp"
53	+	#include "io/StickyAtomTypesSectionParser.hpp"
54	+	#include "io/BondTypesSectionParser.hpp"
55	+	#include "io/BendTypesSectionParser.hpp"
56	+	#include "io/TorsionTypesSectionParser.hpp"
57	+	#include "UseTheForce/ForceFieldCreator.hpp"
58
59	<	#ifdef IS_MPI
17	<	#include "UseTheForce/mpiForceField.h"
18	<	#endif // is_mpi
19	<
20	<
21	<	// define some bond Types
22	<
23	<	#define FIXED_BOND    0
24	<	#define HARMONIC_BOND 1
25	<
26	<
27	<	namespace DUFF_NS {  // restrict the access of the folowing to this file only.
28	<
29	<
30	<	// Declare the structures that will be passed by MPI
31	<
32	<	typedef struct{
33	<	char name[15];
34	<	double mass;
35	<	double epslon;
36	<	double sigma;
37	<	double charge;
38	<	double dipole;
39	<	double w0;
40	<	double v0;
41	<	double v0p;
42	<	double rl;
43	<	double ru;
44	<	double rlp;
45	<	double rup;
46	<	int isSSD;
47	<	int isCharge;
48	<	int isDipole;
49	<	int ident;
50	<	int last;      //  0  -> default
51	<	//  1  -> tells nodes to stop listening
52	<	} atomStruct;
53	<
54	<
55	<	typedef struct{
56	<	char nameA[15];
57	<	char nameB[15];
58	<	double d0;
59	<	double k0;
60	<	int last;      //  0  -> default
61	<	//  1  -> tells nodes to stop listening
62	<	int type;
63	<	} bondStruct;
64	<
65	<
66	<	typedef struct{
67	<	char nameA[15];
68	<	char nameB[15];
69	<	char nameC[15];
70	<	char type[30];
71	<	double k1, k2, k3, t0;
72	<	int last;      //  0  -> default
73	<	//  1  -> tells nodes to stop listening
74	<	} bendStruct;
75	<
76	<
77	<	typedef struct{
78	<	char nameA[15];
79	<	char nameB[15];
80	<	char nameC[15];
81	<	char nameD[15];
82	<	char type[30];
83	<	double k1, k2, k3, k4;
84	<	int last;      //  0  -> default
85	<	//  1  -> tells nodes to stop listening
86	<	} torsionStruct;
87	<
88	<
89	<	int parseAtom(    char *lineBuffer, int lineNum, atomStruct     &info );
90	<	int parseBond(    char *lineBuffer, int lineNum, bondStruct     &info );
91	<	int parseBend(    char *lineBuffer, int lineNum, bendStruct     &info );
92	<	int parseTorsion( char *lineBuffer, int lineNum, torsionStruct  &info );
93	<
94	<
95	<	#ifdef IS_MPI
96	<
97	<	MPI_Datatype mpiAtomStructType;
98	<	MPI_Datatype mpiBondStructType;
99	<	MPI_Datatype mpiBendStructType;
100	<	MPI_Datatype mpiTorsionStructType;
101	<
102	<	#endif
103	<
104	<	class LinkedAtomType {
105	<	public:
106	<	LinkedAtomType(){
107	<	next = NULL;
108	<	name[0] = '\0';
109	<	}
110	<	~LinkedAtomType(){ if( next != NULL ) delete next; }
111	<
112	<	LinkedAtomType* find(char* key){
113	<	if( !strcmp(name, key) ) return this;
114	<	if( next != NULL ) return next->find(key);
115	<	return NULL;
116	<	}
59	>	namespace oopse {
60
61	<	void printMe( void ){
119	<
120	<	std::cerr << "LinkedAtype " << name << ": ident = " << ident << "\n";
121	<	//      if( next != NULL ) next->printMe();
61	>	DUFF::DUFF(){
62
63	<	}
63	>	//set default force field filename
64	>	setForceFieldFileName("DUFF2.frc");
65
66	<	void add( atomStruct &info ){
67	<
68	<	// check for duplicates
69	<
70	<	if( !strcmp( info.name, name ) ){
71	<	sprintf( painCave.errMsg,
72	<	"Duplicate DUFF atom type \"%s\" found in "
73	<	"the DUFF param file./n",
74	<	name );
75	<	painCave.isFatal = 1;
76	<	simError();
77	<	}
78	<
79	<	if( next != NULL ) next->add(info);
80	<	else{
81	<	next = new LinkedAtomType();
82	<	strcpy(next->name, info.name);
83	<	next->isDipole = info.isDipole;
143	<	next->isSSD    = info.isSSD;
144	<	next->mass     = info.mass;
145	<	next->epslon   = info.epslon;
146	<	next->sigma    = info.sigma;
147	<	next->dipole   = info.dipole;
148	<	next->w0       = info.w0;
149	<	next->v0       = info.v0;
150	<	next->v0p      = info.v0p;
151	<	next->rl       = info.rl;
152	<	next->ru       = info.ru;
153	<	next->rlp      = info.rlp;
154	<	next->rup      = info.rup;
155	<	next->ident    = info.ident;
156	<	}
157	<	}
158	<
159	<	#ifdef IS_MPI
66	>	//the order of adding section parsers are important
67	>	//DirectionalAtomTypesSectionParser should be added before AtomTypesSectionParser Since
68	>	//These two section parsers will actually create "real" AtomTypes (AtomTypesSectionParser will create
69	>	//AtomType and DirectionalAtomTypesSectionParser will creat DirectionalAtomType which is a subclass
70	>	//of AtomType, therefore it should come first). Other AtomTypes Section Parser will not create the
71	>	//"real" AtomType, they only add and set some attribute of the AtomType. Thus their order are not
72	>	//important. AtomTypesSectionParser should be added before other atom type section parsers.
73	>	//Make sure they are added after DirectionalAtomTypesSectionParser and AtomTypesSectionParser.
74	>	//The order of BondTypesSectionParser, BendTypesSectionParser and TorsionTypesSectionParser are
75	>	//not important.
76	>	spMan_.push_back(new DirectionalAtomTypesSectionParser());
77	>	spMan_.push_back(new AtomTypesSectionParser());
78	>	spMan_.push_back(new LennardJonesAtomTypesSectionParser());
79	>	spMan_.push_back(new ElectrostaticAtomTypesSectionParser());
80	>	spMan_.push_back(new StickyAtomTypesSectionParser());
81	>	spMan_.push_back(new BondTypesSectionParser());
82	>	spMan_.push_back(new BendTypesSectionParser());
83	>	spMan_.push_back(new TorsionTypesSectionParser());
84
85	<	void duplicate( atomStruct &info ){
162	<	strcpy(info.name, name);
163	<	info.isDipole = isDipole;
164	<	info.isSSD    = isSSD;
165	<	info.mass     = mass;
166	<	info.epslon   = epslon;
167	<	info.sigma    = sigma;
168	<	info.dipole   = dipole;
169	<	info.w0       = w0;
170	<	info.v0       = v0;
171	<	info.v0p      = v0p;
172	<	info.rl       = rl;
173	<	info.ru       = ru;
174	<	info.rlp      = rlp;
175	<	info.rup      = rup;
176	<	info.ident    = ident;
177	<	info.last     = 0;
178	<	}
85	>	}
86
87	+	void DUFF::parse(const std::string& filename) {
88	+	ifstrstream* ffStream;
89	+	ffStream = openForceFieldFile(filename);
90
91	<	#endif
91	>	spMan_.parse(*ffStream, *this);
92
93	<	char name[15];
94	<	int isDipole;
185	<	int isSSD;
186	<	double mass;
187	<	double epslon;
188	<	double sigma;
189	<	double dipole;
190	<	double w0;
191	<	double v0;
192	<	double v0p;
193	<	double rl;
194	<	double ru;
195	<	double rlp;
196	<	double rup;
197	<	int ident;
198	<	LinkedAtomType* next;
199	<	};
93	>	ForceField::AtomTypeContainer::MapTypeIterator i;
94	>	AtomType* at;
95
96	<	class LinkedBondType {
97	<	public:
203	<	LinkedBondType(){
204	<	next = NULL;
205	<	nameA[0] = '\0';
206	<	nameB[0] = '\0';
207	<	}
208	<	~LinkedBondType(){ if( next != NULL ) delete next; }
209	<
210	<	LinkedBondType* find(char* key1, char* key2){
211	<	if( !strcmp(nameA, key1 ) && !strcmp( nameB, key2 ) ) return this;
212	<	if( !strcmp(nameA, key2 ) && !strcmp( nameB, key1 ) ) return this;
213	<	if( next != NULL ) return next->find(key1, key2);
214	<	return NULL;
215	<	}
216	<
217	<
218	<	void add( bondStruct &info ){
219	<
220	<	// check for duplicates
221	<	int dup = 0;
222	<
223	<	if( !strcmp(nameA, info.nameA ) && !strcmp( nameB, info.nameB ) ) dup = 1;
224	<	if( !strcmp(nameA, info.nameB ) && !strcmp( nameB, info.nameA ) ) dup = 1;
225	<
226	<	if(dup){
227	<	sprintf( painCave.errMsg,
228	<	"Duplicate DUFF bond type \"%s - %s\" found in "
229	<	"the DUFF param file./n",
230	<	nameA, nameB );
231	<	painCave.isFatal = 1;
232	<	simError();
233	<	}
234	<
235	<
236	<	if( next != NULL ) next->add(info);
237	<	else{
238	<	next = new LinkedBondType();
239	<	strcpy(next->nameA, info.nameA);
240	<	strcpy(next->nameB, info.nameB);
241	<	next->type = info.type;
242	<	next->d0 = info.d0;
243	<	next->k0 = info.k0;
244	<	}
245	<	}
246	<
247	<	#ifdef IS_MPI
248	<	void duplicate( bondStruct &info ){
249	<	strcpy(info.nameA, nameA);
250	<	strcpy(info.nameB, nameB);
251	<	info.type = type;
252	<	info.d0   = d0;
253	<	info.k0   = k0;
254	<	info.last = 0;
255	<	}
256	<
257	<
258	<	#endif
259	<
260	<	char nameA[15];
261	<	char nameB[15];
262	<	int type;
263	<	double d0;
264	<	double k0;
265	<
266	<	LinkedBondType* next;
267	<	};
268	<
269	<	class LinkedBendType {
270	<	public:
271	<	LinkedBendType(){
272	<	next = NULL;
273	<	nameA[0] = '\0';
274	<	nameB[0] = '\0';
275	<	nameC[0] = '\0';
276	<	type[0] = '\0';
277	<	}
278	<	~LinkedBendType(){ if( next != NULL ) delete next; }
279	<
280	<	LinkedBendType* find( char* key1, char* key2, char* key3 ){
281	<	if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 )
282	<	&& !strcmp( nameC, key3 ) ) return this;
283	<	if( !strcmp( nameA, key3 ) && !strcmp( nameB, key2 )
284	<	&& !strcmp( nameC, key1 ) ) return this;
285	<	if( next != NULL ) return next->find(key1, key2, key3);
286	<	return NULL;
287	<	}
288	<
289	<	void add( bendStruct &info ){
290	<
291	<	// check for duplicates
292	<	int dup = 0;
293	<
294	<	if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB )
295	<	&& !strcmp( nameC, info.nameC ) ) dup = 1;
296	<	if( !strcmp( nameA, info.nameC ) && !strcmp( nameB, info.nameB )
297	<	&& !strcmp( nameC, info.nameA ) ) dup = 1;
298	<
299	<	if(dup){
300	<	sprintf( painCave.errMsg,
301	<	"Duplicate DUFF bend type \"%s - %s - %s\" found in "
302	<	"the DUFF param file./n",
303	<	nameA, nameB, nameC );
304	<	painCave.isFatal = 1;
305	<	simError();
306	<	}
307	<
308	<	if( next != NULL ) next->add(info);
309	<	else{
310	<	next = new LinkedBendType();
311	<	strcpy(next->nameA, info.nameA);
312	<	strcpy(next->nameB, info.nameB);
313	<	strcpy(next->nameC, info.nameC);
314	<	strcpy(next->type,  info.type);
315	<	next->k1 = info.k1;
316	<	next->k2 = info.k2;
317	<	next->k3 = info.k3;
318	<	next->t0 = info.t0;
319	<	}
320	<	}
321	<
322	<	#ifdef IS_MPI
323	<
324	<	void duplicate( bendStruct &info ){
325	<	strcpy(info.nameA, nameA);
326	<	strcpy(info.nameB, nameB);
327	<	strcpy(info.nameC, nameC);
328	<	strcpy(info.type,  type);
329	<	info.k1   = k1;
330	<	info.k2   = k2;
331	<	info.k3   = k3;
332	<	info.t0   = t0;
333	<	info.last = 0;
334	<	}
335	<
336	<	#endif // is_mpi
337	<
338	<	char nameA[15];
339	<	char nameB[15];
340	<	char nameC[15];
341	<	char type[30];
342	<	double k1, k2, k3, t0;
343	<
344	<	LinkedBendType* next;
345	<	};
346	<
347	<	class LinkedTorsionType {
348	<	public:
349	<	LinkedTorsionType(){
350	<	next = NULL;
351	<	nameA[0] = '\0';
352	<	nameB[0] = '\0';
353	<	nameC[0] = '\0';
354	<	type[0] = '\0';
355	<	}
356	<	~LinkedTorsionType(){ if( next != NULL ) delete next; }
357	<
358	<	LinkedTorsionType* find( char* key1, char* key2, char* key3, char* key4 ){
359	<
360	<
361	<
362	<
363	<	if( !strcmp( nameA, key1 ) && !strcmp( nameB, key2 ) &&
364	<	!strcmp( nameC, key3 ) && !strcmp( nameD, key4 ) ) return this;
365	<
366	<	if( !strcmp( nameA, key4 ) && !strcmp( nameB, key3 ) &&
367	<	!strcmp( nameC, key2 ) && !strcmp( nameD, key1 ) ) return this;
368	<
369	<	if( next != NULL ) return next->find(key1, key2, key3, key4);
370	<	return NULL;
371	<	}
372	<
373	<	void add( torsionStruct &info ){
374	<
375	<	// check for duplicates
376	<	int dup = 0;
377	<
378	<	if( !strcmp( nameA, info.nameA ) && !strcmp( nameB, info.nameB ) &&
379	<	!strcmp( nameC, info.nameC ) && !strcmp( nameD, info.nameD ) ) dup = 1;
380	<
381	<	if( !strcmp( nameA, info.nameD ) && !strcmp( nameB, info.nameC ) &&
382	<	!strcmp( nameC, info.nameB ) && !strcmp( nameD, info.nameA ) ) dup = 1;
383	<
384	<	if(dup){
385	<	sprintf( painCave.errMsg,
386	<	"Duplicate DUFF torsion type \"%s - %s - %s - %s\" found in "
387	<	"the DUFF param file./n", nameA, nameB, nameC, nameD );
388	<	painCave.isFatal = 1;
389	<	simError();
390	<	}
391	<
392	<	if( next != NULL ) next->add(info);
393	<	else{
394	<	next = new LinkedTorsionType();
395	<	strcpy(next->nameA, info.nameA);
396	<	strcpy(next->nameB, info.nameB);
397	<	strcpy(next->nameC, info.nameC);
398	<	strcpy(next->nameD, info.nameD);
399	<	strcpy(next->type,  info.type);
400	<	next->k1 = info.k1;
401	<	next->k2 = info.k2;
402	<	next->k3 = info.k3;
403	<	next->k4 = info.k4;
404	<
405	<	}
406	<	}
407	<
408	<	#ifdef IS_MPI
409	<
410	<	void duplicate( torsionStruct &info ){
411	<	strcpy(info.nameA, nameA);
412	<	strcpy(info.nameB, nameB);
413	<	strcpy(info.nameC, nameC);
414	<	strcpy(info.nameD, nameD);
415	<	strcpy(info.type,  type);
416	<	info.k1   = k1;
417	<	info.k2   = k2;
418	<	info.k3   = k3;
419	<	info.k4   = k4;
420	<	info.last = 0;
96	>	for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
97	>	at->makeFortranAtomType();
98		}
99
100	<	#endif
101	<
425	<	char nameA[15];
426	<	char nameB[15];
427	<	char nameC[15];
428	<	char nameD[15];
429	<	char type[30];
430	<	double k1, k2, k3, k4;
431	<
432	<	LinkedTorsionType* next;
433	<	};
434	<
435	<
436	<	LinkedAtomType* headAtomType;
437	<	LinkedAtomType* currentAtomType;
438	<	LinkedBondType* headBondType;
439	<	LinkedBondType* currentBondType;
440	<	LinkedBendType* headBendType;
441	<	LinkedBendType* currentBendType;
442	<	LinkedTorsionType* headTorsionType;
443	<	LinkedTorsionType* currentTorsionType;
444	<
445	<	} // namespace
446	<
447	<	using namespace DUFF_NS;
448	<
449	<
450	<	//****************************************************************
451	<	// begins the actual forcefield stuff.
452	<	//****************************************************************
453	<
454	<
455	<	DUFF::DUFF(){
456	<
457	<	char fileName[200];
458	<	char* ffPath_env = "FORCE_PARAM_PATH";
459	<	char* ffPath;
460	<	char temp[200];
461	<
462	<	headAtomType       = NULL;
463	<	currentAtomType    = NULL;
464	<	headBondType       = NULL;
465	<	currentBondType    = NULL;
466	<	headBendType       = NULL;
467	<	currentBendType    = NULL;
468	<	headTorsionType    = NULL;
469	<	currentTorsionType = NULL;
470	<
471	<
472	<	#ifdef IS_MPI
473	<	int i;
474	<
475	<	// **********************************************************************
476	<	// Init the atomStruct mpi type
477	<
478	<	atomStruct atomProto; // mpiPrototype
479	<	int atomBC[3] = {15,12,5};  // block counts
480	<	MPI_Aint atomDspls[3];           // displacements
481	<	MPI_Datatype atomMbrTypes[3];    // member mpi types
482	<
483	<	MPI_Address(&atomProto.name, &atomDspls[0]);
484	<	MPI_Address(&atomProto.mass, &atomDspls[1]);
485	<	MPI_Address(&atomProto.isSSD, &atomDspls[2]);
486	<
487	<	atomMbrTypes[0] = MPI_CHAR;
488	<	atomMbrTypes[1] = MPI_DOUBLE;
489	<	atomMbrTypes[2] = MPI_INT;
490	<
491	<	for (i=2; i >= 0; i--) atomDspls[i] -= atomDspls[0];
492	<
493	<	MPI_Type_struct(3, atomBC, atomDspls, atomMbrTypes, &mpiAtomStructType);
494	<	MPI_Type_commit(&mpiAtomStructType);
495	<
496	<
497	<	// **********************************************************************
498	<	// Init the bondStruct mpi type
499	<
500	<	bondStruct bondProto; // mpiPrototype
501	<	int bondBC[3] = {30,2,2};  // block counts
502	<	MPI_Aint bondDspls[3];           // displacements
503	<	MPI_Datatype bondMbrTypes[3];    // member mpi types
504	<
505	<	MPI_Address(&bondProto.nameA, &bondDspls[0]);
506	<	MPI_Address(&bondProto.d0,    &bondDspls[1]);
507	<	MPI_Address(&bondProto.last,  &bondDspls[2]);
508	<
509	<	bondMbrTypes[0] = MPI_CHAR;
510	<	bondMbrTypes[1] = MPI_DOUBLE;
511	<	bondMbrTypes[2] = MPI_INT;
512	<
513	<	for (i=2; i >= 0; i--) bondDspls[i] -= bondDspls[0];
514	<
515	<	MPI_Type_struct(3, bondBC, bondDspls, bondMbrTypes, &mpiBondStructType);
516	<	MPI_Type_commit(&mpiBondStructType);
517	<
518	<
519	<	// **********************************************************************
520	<	// Init the bendStruct mpi type
521	<
522	<	bendStruct bendProto; // mpiPrototype
523	<	int bendBC[3] = {75,4,1};  // block counts
524	<	MPI_Aint bendDspls[3];           // displacements
525	<	MPI_Datatype bendMbrTypes[3];    // member mpi types
526	<
527	<	MPI_Address(&bendProto.nameA, &bendDspls[0]);
528	<	MPI_Address(&bendProto.k1,    &bendDspls[1]);
529	<	MPI_Address(&bendProto.last,  &bendDspls[2]);
530	<
531	<	bendMbrTypes[0] = MPI_CHAR;
532	<	bendMbrTypes[1] = MPI_DOUBLE;
533	<	bendMbrTypes[2] = MPI_INT;
534	<
535	<	for (i=2; i >= 0; i--) bendDspls[i] -= bendDspls[0];
536	<
537	<	MPI_Type_struct(3, bendBC, bendDspls, bendMbrTypes, &mpiBendStructType);
538	<	MPI_Type_commit(&mpiBendStructType);
539	<
540	<
541	<	// **********************************************************************
542	<	// Init the torsionStruct mpi type
543	<
544	<	torsionStruct torsionProto; // mpiPrototype
545	<	int torsionBC[3] = {90,4,1};  // block counts
546	<	MPI_Aint torsionDspls[3];           // displacements
547	<	MPI_Datatype torsionMbrTypes[3];    // member mpi types
548	<
549	<	MPI_Address(&torsionProto.nameA, &torsionDspls[0]);
550	<	MPI_Address(&torsionProto.k1,    &torsionDspls[1]);
551	<	MPI_Address(&torsionProto.last,  &torsionDspls[2]);
552	<
553	<	torsionMbrTypes[0] = MPI_CHAR;
554	<	torsionMbrTypes[1] = MPI_DOUBLE;
555	<	torsionMbrTypes[2] = MPI_INT;
556	<
557	<	for (i=2; i >= 0; i--) torsionDspls[i] -= torsionDspls[0];
558	<
559	<	MPI_Type_struct(3, torsionBC, torsionDspls, torsionMbrTypes,
560	<	&mpiTorsionStructType);
561	<	MPI_Type_commit(&mpiTorsionStructType);
562	<
563	<	// ***********************************************************************
564	<
565	<	if( worldRank == 0 ){
566	<	#endif
567	<
568	<	// generate the force file name
569	<
570	<	strcpy( fileName, "DUFF.frc" );
571	<	//    fprintf( stderr,"Trying to open %s\n", fileName );
572	<
573	<	// attempt to open the file in the current directory first.
574	<
575	<	frcFile = fopen( fileName, "r" );
576	<
577	<	if( frcFile == NULL ){
578	<
579	<	// next see if the force path enviorment variable is set
580	<
581	<	ffPath = getenv( ffPath_env );
582	<	if( ffPath == NULL ) {
583	<	STR_DEFINE(ffPath, FRC_PATH );
584	<	}
585	<
586	<
587	<	strcpy( temp, ffPath );
588	<	strcat( temp, "/" );
589	<	strcat( temp, fileName );
590	<	strcpy( fileName, temp );
591	<
592	<	frcFile = fopen( fileName, "r" );
593	<
594	<	if( frcFile == NULL ){
595	<
596	<	sprintf( painCave.errMsg,
597	<	"Error opening the force field parameter file:\n"
598	<	"\t%s\n"
599	<	"\tHave you tried setting the FORCE_PARAM_PATH environment "
600	<	"variable?\n",
601	<	fileName );
602	<	painCave.severity = OOPSE_ERROR;
603	<	painCave.isFatal = 1;
604	<	simError();
605	<	}
100	>	for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
101	>	at->complete();
102		}
103
608	–	#ifdef IS_MPI
609	–	}
610	–
611	–	sprintf( checkPointMsg, "DUFF file opened sucessfully." );
612	–	MPIcheckPoint();
613	–
614	–	#endif // is_mpi
615	–	}
616	–
617	–
618	–	DUFF::~DUFF(){
619	–
620	–	if( headAtomType != NULL ) delete headAtomType;
621	–	if( headBondType != NULL ) delete headBondType;
622	–	if( headBendType != NULL ) delete headBendType;
623	–	if( headTorsionType != NULL ) delete headTorsionType;
624	–
625	–	#ifdef IS_MPI
626	–	if( worldRank == 0 ){
627	–	#endif // is_mpi
628	–
629	–	fclose( frcFile );
630	–
631	–	#ifdef IS_MPI
632	–	}
633	–	#endif // is_mpi
634	–	}
635	–
636	–	void DUFF::cleanMe( void ){
637	–
638	–	#ifdef IS_MPI
639	–
640	–	// keep the linked lists in the mpi version
641	–
642	–	#else // is_mpi
643	–
644	–	// delete the linked lists in the single processor version
645	–
646	–	if( headAtomType != NULL ) delete headAtomType;
647	–	if( headBondType != NULL ) delete headBondType;
648	–	if( headBendType != NULL ) delete headBendType;
649	–	if( headTorsionType != NULL ) delete headTorsionType;
650	–
651	–	#endif // is_mpi
652	–	}
653	–
654	–
655	–	void DUFF::initForceField(){
656	–
657	–	initFortran( entry_plug->useReactionField );
658	–	}
659	–
660	–
661	–	void DUFF::readParams( void ){
662	–
663	–	int identNum, isError;
664	–
665	–	atomStruct atomInfo;
666	–	bondStruct bondInfo;
667	–	bendStruct bendInfo;
668	–	torsionStruct torsionInfo;
669	–
670	–	AtomType* at;
671	–
672	–	bigSigma = 0.0;
673	–
674	–	atomInfo.last = 1;
675	–	bondInfo.last = 1;
676	–	bendInfo.last = 1;
677	–	torsionInfo.last = 1;
678	–
679	–	// read in the atom info
680	–
681	–	#ifdef IS_MPI
682	–	if( worldRank == 0 ){
683	–	#endif
684	–
685	–	// read in the atom types.
686	–
687	–	headAtomType = new LinkedAtomType;
688	–
689	–	fastForward( "AtomTypes", "initializeAtoms" );
690	–
691	–	// we are now at the AtomTypes section.
692	–
693	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
694	–	lineNum++;
695	–
696	–
697	–	// read a line, and start parseing out the atom types
698	–
699	–	if( eof_test == NULL ){
700	–	sprintf( painCave.errMsg,
701	–	"Error in reading Atoms from force file at line %d.\n",
702	–	lineNum );
703	–	painCave.isFatal = 1;
704	–	simError();
705	–	}
706	–
707	–	identNum = 1;
708	–	// stop reading at end of file, or at next section
709	–	while( readLine[0] != '#' && eof_test != NULL ){
710	–
711	–	// toss comment lines
712	–	if( readLine[0] != '!' ){
713	–
714	–	// the parser returns 0 if the line was blank
715	–	if( parseAtom( readLine, lineNum, atomInfo ) ){
716	–	atomInfo.ident = identNum;
717	–	headAtomType->add( atomInfo );;
718	–	identNum++;
719	–	}
720	–	}
721	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
722	–	lineNum++;
723	–	}
724	–
725	–	#ifdef IS_MPI
726	–
727	–	// send out the linked list to all the other processes
728	–
729	–	sprintf( checkPointMsg,
730	–	"DUFF atom structures read successfully." );
731	–	MPIcheckPoint();
732	–
733	–	currentAtomType = headAtomType->next; //skip the first element who is a place holder.
734	–	while( currentAtomType != NULL ){
735	–	currentAtomType->duplicate( atomInfo );
736	–
737	–	sendFrcStruct( &atomInfo, mpiAtomStructType );
738	–
739	–	sprintf( checkPointMsg,
740	–	"successfully sent DUFF force type: \"%s\"\n",
741	–	atomInfo.name );
742	–	MPIcheckPoint();
743	–
744	–	currentAtomType = currentAtomType->next;
745	–	}
746	–	atomInfo.last = 1;
747	–	sendFrcStruct( &atomInfo, mpiAtomStructType );
748	–
749	–	}
750	–
751	–	else{
752	–
753	–	// listen for node 0 to send out the force params
754	–
755	–	MPIcheckPoint();
756	–
757	–	headAtomType = new LinkedAtomType;
758	–	receiveFrcStruct( &atomInfo, mpiAtomStructType );
759	–
760	–	while( !atomInfo.last ){
761	–
762	–	headAtomType->add( atomInfo );
763	–
764	–	MPIcheckPoint();
765	–
766	–	receiveFrcStruct( &atomInfo, mpiAtomStructType );
767	–	}
768	–	}
769	–
770	–	#endif // is_mpi
771	–
772	–	// dummy variables
773	–
774	–	currentAtomType = headAtomType->next;;
775	–	while( currentAtomType != NULL ){
776	–
777	–	if( currentAtomType->name[0] != '\0' ){
778	–
779	–	if (currentAtomType->isSSD || currentAtomType->isDipole)
780	–	DirectionalAtomType* at = new DirectionalAtomType();
781	–	else
782	–	AtomType* at = new AtomType();
783	–
784	–	if (currentAtomType->isSSD) {
785	–	((DirectionalAtomType*)at)->setSticky();
786	–	entry_plug->useSticky = 1;
787	–	}
788	–
789	–	if (currentAtomType->isDipole) {
790	–	((DirectionalAtomType*)at)->setDipole();
791	–	entry_plug->useDipoles = 1;
792	–	}
793	–
794	–	at->setIdent(currentAtomType->ident);
795	–	at->setName(currentAtomType->name);
796	–	at->setLennardJones();
797	–	at->complete();
798	–	}
799	–	currentAtomType = currentAtomType->next;
800	–	}
801	–
802	–	currentAtomType = headAtomType->next;;
803	–	while( currentAtomType != NULL ){
804	–
805	–	if( currentAtomType->name[0] != '\0' ){
806	–	isError = 0;
807	–	newLJtype( &(currentAtomType->ident), &(currentAtomType->sigma),
808	–	&(currentAtomType->epslon), &isError);
809	–	if( isError ){
810	–	sprintf( painCave.errMsg,
811	–	"Error initializing the \"%s\" LJ type in fortran\n",
812	–	currentAtomType->name );
813	–	painCave.isFatal = 1;
814	–	simError();
815	–	}
816	–
817	–	if (currentAtomType->isDipole) {
818	–	newDipoleType(&(currentAtomType->ident), &(currentAtomType->dipole),
819	–	&isError);
820	–	if( isError ){
821	–	sprintf( painCave.errMsg,
822	–	"Error initializing the \"%s\" dipole type in fortran\n",
823	–	currentAtomType->name );
824	–	painCave.isFatal = 1;
825	–	simError();
826	–	}
827	–	}
828	–
829	–	if(currentAtomType->isSSD) {
830	–	makeStickyType( &(currentAtomType->w0), &(currentAtomType->v0),
831	–	&(currentAtomType->v0p),
832	–	&(currentAtomType->rl), &(currentAtomType->ru),
833	–	&(currentAtomType->rlp), &(currentAtomType->rup));
834	–	}
835	–
836	–	}
837	–	currentAtomType = currentAtomType->next;
838	–	}
839	–
840	–	#ifdef IS_MPI
841	–	sprintf( checkPointMsg,
842	–	"DUFF atom structures successfully sent to fortran\n" );
843	–	MPIcheckPoint();
844	–	#endif // is_mpi
845	–
846	–
847	–
848	–	// read in the bonds
849	–
850	–	#ifdef IS_MPI
851	–	if( worldRank == 0 ){
852	–	#endif
853	–
854	–	// read in the bond types.
855	–
856	–	headBondType = new LinkedBondType;
857	–
858	–	fastForward( "BondTypes", "initializeBonds" );
859	–
860	–	// we are now at the bondTypes section
861	–
862	–	eof_test =  fgets( readLine, sizeof(readLine), frcFile );
863	–	lineNum++;
864	–
865	–
866	–	// read a line, and start parseing out the atom types
867	–
868	–	if( eof_test == NULL ){
869	–	sprintf( painCave.errMsg,
870	–	"Error in reading bonds from force file at line %d.\n",
871	–	lineNum );
872	–	painCave.isFatal = 1;
873	–	simError();
874	–	}
875	–
876	–	// stop reading at end of file, or at next section
877	–	while( readLine[0] != '#' && eof_test != NULL ){
878	–
879	–	// toss comment lines
880	–	if( readLine[0] != '!' ){
881	–
882	–	// the parser returns 0 if the line was blank
883	–	if( parseBond( readLine, lineNum, bondInfo ) ){
884	–	headBondType->add( bondInfo );
885	–	}
886	–	}
887	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
888	–	lineNum++;
889	–	}
890	–
891	–	#ifdef IS_MPI
892	–
893	–	// send out the linked list to all the other processes
894	–
895	–	sprintf( checkPointMsg,
896	–	"DUFF bond structures read successfully." );
897	–	MPIcheckPoint();
898	–
899	–	currentBondType = headBondType->next;
900	–	while( currentBondType != NULL ){
901	–	currentBondType->duplicate( bondInfo );
902	–	sendFrcStruct( &bondInfo, mpiBondStructType );
903	–	currentBondType = currentBondType->next;
904	–	}
905	–	bondInfo.last = 1;
906	–	sendFrcStruct( &bondInfo, mpiBondStructType );
907	–
908	–	}
909	–
910	–	else{
911	–
912	–	// listen for node 0 to send out the force params
913	–
914	–	MPIcheckPoint();
915	–
916	–	headBondType = new LinkedBondType;
917	–	receiveFrcStruct( &bondInfo, mpiBondStructType );
918	–	while( !bondInfo.last ){
919	–
920	–	headBondType->add( bondInfo );
921	–	receiveFrcStruct( &bondInfo, mpiBondStructType );
922	–	}
923	–	}
924	–
925	–	sprintf( checkPointMsg,
926	–	"DUFF bond structures broadcast successfully." );
927	–	MPIcheckPoint();
928	–
929	–	#endif // is_mpi
930	–
931	–
932	–	// read in the bends
933	–
934	–	#ifdef IS_MPI
935	–	if( worldRank == 0 ){
936	–	#endif
937	–
938	–	// read in the bend types.
939	–
940	–	headBendType = new LinkedBendType;
941	–
942	–	fastForward( "BendTypes", "initializeBends" );
943	–
944	–	// we are now at the bendTypes section
945	–
946	–	eof_test =  fgets( readLine, sizeof(readLine), frcFile );
947	–	lineNum++;
948	–
949	–	// read a line, and start parseing out the bend types
950	–
951	–	if( eof_test == NULL ){
952	–	sprintf( painCave.errMsg,
953	–	"Error in reading bends from force file at line %d.\n",
954	–	lineNum );
955	–	painCave.isFatal = 1;
956	–	simError();
957	–	}
958	–
959	–	// stop reading at end of file, or at next section
960	–	while( readLine[0] != '#' && eof_test != NULL ){
961	–
962	–	// toss comment lines
963	–	if( readLine[0] != '!' ){
964	–
965	–	// the parser returns 0 if the line was blank
966	–	if( parseBend( readLine, lineNum, bendInfo ) ){
967	–	headBendType->add( bendInfo );
968	–	}
969	–	}
970	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
971	–	lineNum++;
972	–	}
973	–
974	–	#ifdef IS_MPI
975	–
976	–	// send out the linked list to all the other processes
977	–
978	–	sprintf( checkPointMsg,
979	–	"DUFF bend structures read successfully." );
980	–	MPIcheckPoint();
981	–
982	–	currentBendType = headBendType->next;
983	–	while( currentBendType != NULL ){
984	–	currentBendType->duplicate( bendInfo );
985	–	sendFrcStruct( &bendInfo, mpiBendStructType );
986	–	currentBendType = currentBendType->next;
987	–	}
988	–	bendInfo.last = 1;
989	–	sendFrcStruct( &bendInfo, mpiBendStructType );
990	–
991	–	}
992	–
993	–	else{
994	–
995	–	// listen for node 0 to send out the force params
996	–
997	–	MPIcheckPoint();
998	–
999	–	headBendType = new LinkedBendType;
1000	–	receiveFrcStruct( &bendInfo, mpiBendStructType );
1001	–	while( !bendInfo.last ){
1002	–
1003	–	headBendType->add( bendInfo );
1004	–	receiveFrcStruct( &bendInfo, mpiBendStructType );
1005	–	}
1006	–	}
1007	–
1008	–	sprintf( checkPointMsg,
1009	–	"DUFF bend structures broadcast successfully." );
1010	–	MPIcheckPoint();
1011	–
1012	–	#endif // is_mpi
1013	–
1014	–
1015	–	// read in the torsions
1016	–
1017	–	#ifdef IS_MPI
1018	–	if( worldRank == 0 ){
1019	–	#endif
1020	–
1021	–	// read in the torsion types.
1022	–
1023	–	headTorsionType = new LinkedTorsionType;
1024	–
1025	–	fastForward( "TorsionTypes", "initializeTorsions" );
1026	–
1027	–	// we are now at the torsionTypes section
1028	–
1029	–	eof_test =  fgets( readLine, sizeof(readLine), frcFile );
1030	–	lineNum++;
1031	–
1032	–
1033	–	// read a line, and start parseing out the atom types
1034	–
1035	–	if( eof_test == NULL ){
1036	–	sprintf( painCave.errMsg,
1037	–	"Error in reading torsions from force file at line %d.\n",
1038	–	lineNum );
1039	–	painCave.isFatal = 1;
1040	–	simError();
1041	–	}
1042	–
1043	–	// stop reading at end of file, or at next section
1044	–	while( readLine[0] != '#' && eof_test != NULL ){
1045	–
1046	–	// toss comment lines
1047	–	if( readLine[0] != '!' ){
1048	–
1049	–	// the parser returns 0 if the line was blank
1050	–	if( parseTorsion( readLine, lineNum, torsionInfo ) ){
1051	–	headTorsionType->add( torsionInfo );
1052	–
1053	–	}
1054	–	}
1055	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1056	–	lineNum++;
1057	–	}
1058	–
1059	–	#ifdef IS_MPI
1060	–
1061	–	// send out the linked list to all the other processes
1062	–
1063	–	sprintf( checkPointMsg,
1064	–	"DUFF torsion structures read successfully." );
1065	–	MPIcheckPoint();
1066	–
1067	–	currentTorsionType = headTorsionType->next;
1068	–	while( currentTorsionType != NULL ){
1069	–	currentTorsionType->duplicate( torsionInfo );
1070	–	sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1071	–	currentTorsionType = currentTorsionType->next;
1072	–	}
1073	–	torsionInfo.last = 1;
1074	–	sendFrcStruct( &torsionInfo, mpiTorsionStructType );
1075	–
1076	–	}
1077	–
1078	–	else{
1079	–
1080	–	// listen for node 0 to send out the force params
1081	–
1082	–	MPIcheckPoint();
1083	–
1084	–	headTorsionType = new LinkedTorsionType;
1085	–	receiveFrcStruct( &torsionInfo, mpiTorsionStructType );
1086	–	while( !torsionInfo.last ){
1087	–
1088	–	headTorsionType->add( torsionInfo );
1089	–	receiveFrcStruct( &torsionInfo, mpiTorsionStructType );
1090	–	}
1091	–	}
1092	–
1093	–	sprintf( checkPointMsg,
1094	–	"DUFF torsion structures broadcast successfully." );
1095	–	MPIcheckPoint();
1096	–
1097	–	#endif // is_mpi
1098	–
1099	–	entry_plug->useLennardJones = 1;
1100	–	}
1101	–
1102	–
1103	–
1104	–	void DUFF::initializeAtoms( int nAtoms, Atom** the_atoms ){
1105	–
1106	–
1107	–	//////////////////////////////////////////////////
1108	–	// a quick water fix
1109	–
1110	–	double waterI[3][3];
1111	–	waterI[0][0] = 1.76958347772500;
1112	–	waterI[0][1] = 0.0;
1113	–	waterI[0][2] = 0.0;
1114	–
1115	–	waterI[1][0] = 0.0;
1116	–	waterI[1][1] = 0.614537057924513;
1117	–	waterI[1][2] = 0.0;
1118	–
1119	–	waterI[2][0] = 0.0;
1120	–	waterI[2][1] = 0.0;
1121	–	waterI[2][2] = 1.15504641980049;
1122	–
1123	–
1124	–	double headI[3][3];
1125	–	headI[0][0] = 1125;
1126	–	headI[0][1] = 0.0;
1127	–	headI[0][2] = 0.0;
1128	–
1129	–	headI[1][0] = 0.0;
1130	–	headI[1][1] = 1125;
1131	–	headI[1][2] = 0.0;
1132	–
1133	–	headI[2][0] = 0.0;
1134	–	headI[2][1] = 0.0;
1135	–	headI[2][2] = 250;
1136	–
1137	–	//////////////////////////////////////////////////
1138	–
1139	–
1140	–	// initialize the atoms
1141	–
1142	–	DirectionalAtom* dAtom;
1143	–	double ji[3];
1144	–
1145	–	for(int i=0; i<nAtoms; i++ ){
1146	–
1147	–	currentAtomType = headAtomType->find( the_atoms[i]->getType() );
1148	–	if( currentAtomType == NULL ){
1149	–	sprintf( painCave.errMsg,
1150	–	"AtomType error, %s not found in force file.\n",
1151	–	the_atoms[i]->getType() );
1152	–	painCave.isFatal = 1;
1153	–	simError();
1154	–	}
1155	–
1156	–	the_atoms[i]->setMass( currentAtomType->mass );
1157	–	the_atoms[i]->setIdent( currentAtomType->ident );
1158	–
1159	–	if( bigSigma < currentAtomType->sigma ) bigSigma = currentAtomType->sigma;
1160	–
1161	–	if( currentAtomType->isDipole ){
1162	–	if( the_atoms[i]->isDirectional() ){
1163	–
1164	–	dAtom = (DirectionalAtom *) the_atoms[i];
1165	–	dAtom->setHasDipole( 1 );
1166	–
1167	–	ji[0] = 0.0;
1168	–	ji[1] = 0.0;
1169	–	ji[2] = 0.0;
1170	–
1171	–	dAtom->setJ( ji );
1172	–
1173	–	if(!strcmp("SSD",the_atoms[i]->getType())){
1174	–	dAtom->setI( waterI );
1175	–	}
1176	–	else if(!strcmp("HEAD",the_atoms[i]->getType())){
1177	–	dAtom->setI( headI );
1178	–	}
1179	–	else{
1180	–	sprintf(painCave.errMsg,
1181	–	"AtmType error, %s does not have a moment of inertia set.\n",
1182	–	the_atoms[i]->getType() );
1183	–	painCave.isFatal = 1;
1184	–	simError();
1185	–	}
1186	–	entry_plug->n_dipoles++;
1187	–	}
1188	–	else{
1189	–
1190	–	sprintf( painCave.errMsg,
1191	–	"DUFF error: Atom \"%s\" is a dipole, yet no standard"
1192	–	" orientation was specifed in the meta-data file.\n",
1193	–	currentAtomType->name );
1194	–	painCave.isFatal = 1;
1195	–	simError();
1196	–	}
1197	–	}
1198	–	else{
1199	–	if( the_atoms[i]->isDirectional() ){
1200	–	sprintf( painCave.errMsg,
1201	–	"DUFF error: Atom \"%s\" was given a standard "
1202	–	"orientation in the meta-data file, yet it is not a dipole.\n",
1203	–	currentAtomType->name);
1204	–	painCave.isFatal = 1;
1205	–	simError();
1206	–	}
1207	–	}
1208	–	}
1209	–	}
1210	–
1211	–	void DUFF::initializeBonds( int nBonds, Bond** bondArray,
1212	–	bond_pair* the_bonds ){
1213	–	int i,a,b;
1214	–	char* atomA;
1215	–	char* atomB;
1216	–
1217	–	Atom** the_atoms;
1218	–	the_atoms = entry_plug->atoms;
1219	–
1220	–
1221	–	// initialize the Bonds
1222	–
1223	–	for( i=0; i<nBonds; i++ ){
1224	–
1225	–	a = the_bonds[i].a;
1226	–	b = the_bonds[i].b;
1227	–
1228	–	atomA = the_atoms[a]->getType();
1229	–	atomB = the_atoms[b]->getType();
1230	–	currentBondType = headBondType->find( atomA, atomB );
1231	–	if( currentBondType == NULL ){
1232	–	sprintf( painCave.errMsg,
1233	–	"BondType error, %s - %s not found in force file.\n",
1234	–	atomA, atomB );
1235	–	painCave.isFatal = 1;
1236	–	simError();
1237	–	}
1238	–
1239	–	switch( currentBondType->type ){
1240	–
1241	–	case FIXED_BOND:
1242	–
1243	–	bondArray[i] = new ConstrainedBond( *the_atoms[a],
1244	–	*the_atoms[b],
1245	–	currentBondType->d0 );
1246	–	entry_plug->n_constraints++;
1247	–	break;
1248	–
1249	–	case HARMONIC_BOND:
1250	–
1251	–	bondArray[i] = new HarmonicBond( *the_atoms[a],
1252	–	*the_atoms[b],
1253	–	currentBondType->d0,
1254	–	currentBondType->k0 );
1255	–	break;
1256	–
1257	–	default:
1258	–
1259	–	break;
1260	–	// do nothing
1261	–	}
1262	–	}
1263	–	}
1264	–
1265	–	void DUFF::initializeBends( int nBends, Bend** bendArray,
1266	–	bend_set* the_bends ){
1267	–
1268	–	QuadraticBend* qBend;
1269	–	GhostBend* gBend;
1270	–	Atom** the_atoms;
1271	–	the_atoms = entry_plug->atoms;
1272	–
1273	–	int i, a, b, c;
1274	–	char* atomA;
1275	–	char* atomB;
1276	–	char* atomC;
1277	–
1278	–	// initialize the Bends
1279	–
1280	–	for( i=0; i<nBends; i++ ){
1281	–
1282	–	a = the_bends[i].a;
1283	–	b = the_bends[i].b;
1284	–	c = the_bends[i].c;
1285	–
1286	–	atomA = the_atoms[a]->getType();
1287	–	atomB = the_atoms[b]->getType();
1288	–
1289	–	if( the_bends[i].isGhost ) atomC = "GHOST";
1290	–	else atomC = the_atoms[c]->getType();
1291	–
1292	–	currentBendType = headBendType->find( atomA, atomB, atomC );
1293	–	if( currentBendType == NULL ){
1294	–	sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
1295	–	" in force file.\n",
1296	–	atomA, atomB, atomC );
1297	–	painCave.isFatal = 1;
1298	–	simError();
1299	–	}
1300	–
1301	–	if( !strcmp( currentBendType->type, "quadratic" ) ){
1302	–
1303	–	if( the_bends[i].isGhost){
1304	–
1305	–	if( the_bends[i].ghost == b ){
1306	–	// do nothing
1307	–	}
1308	–	else if( the_bends[i].ghost == a ){
1309	–	c = a;
1310	–	a = b;
1311	–	b = c;
1312	–	}
1313	–	else{
1314	–	sprintf( painCave.errMsg,
1315	–	"BendType error, %s - %s - %s,\n"
1316	–	"  --> central atom is not "
1317	–	"correctly identified with the "
1318	–	"\"ghostVectorSource = \" tag.\n",
1319	–	atomA, atomB, atomC );
1320	–	painCave.isFatal = 1;
1321	–	simError();
1322	–	}
1323	–
1324	–	gBend = new GhostBend( *the_atoms[a],
1325	–	*the_atoms[b]);
1326	–
1327	–	gBend->setConstants( currentBendType->k1,
1328	–	currentBendType->k2,
1329	–	currentBendType->k3,
1330	–	currentBendType->t0 );
1331	–	bendArray[i] = gBend;
1332	–	}
1333	–	else{
1334	–	qBend = new QuadraticBend( *the_atoms[a],
1335	–	*the_atoms[b],
1336	–	*the_atoms[c] );
1337	–	qBend->setConstants( currentBendType->k1,
1338	–	currentBendType->k2,
1339	–	currentBendType->k3,
1340	–	currentBendType->t0 );
1341	–	bendArray[i] = qBend;
1342	–	}
1343	–	}
1344	–	}
1345	–	}
1346	–
1347	–	void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1348	–	torsion_set* the_torsions ){
1349	–
1350	–	int i, a, b, c, d;
1351	–	char* atomA;
1352	–	char* atomB;
1353	–	char* atomC;
1354	–	char* atomD;
1355	–
1356	–	CubicTorsion* cTors;
1357	–	Atom** the_atoms;
1358	–	the_atoms = entry_plug->atoms;
1359	–
1360	–	// initialize the Torsions
1361	–
1362	–	for( i=0; i<nTorsions; i++ ){
1363	–
1364	–	a = the_torsions[i].a;
1365	–	b = the_torsions[i].b;
1366	–	c = the_torsions[i].c;
1367	–	d = the_torsions[i].d;
1368	–
1369	–	atomA = the_atoms[a]->getType();
1370	–	atomB = the_atoms[b]->getType();
1371	–	atomC = the_atoms[c]->getType();
1372	–	atomD = the_atoms[d]->getType();
1373	–	currentTorsionType = headTorsionType->find( atomA, atomB, atomC, atomD );
1374	–	if( currentTorsionType == NULL ){
1375	–	sprintf( painCave.errMsg,
1376	–	"TorsionType error, %s - %s - %s - %s not found"
1377	–	" in force file.\n",
1378	–	atomA, atomB, atomC, atomD );
1379	–	painCave.isFatal = 1;
1380	–	simError();
1381	–	}
1382	–
1383	–	if( !strcmp( currentTorsionType->type, "cubic" ) ){
1384	–
1385	–	cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1386	–	*the_atoms[c], *the_atoms[d] );
1387	–	cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1388	–	currentTorsionType->k3, currentTorsionType->k4 );
1389	–	torsionArray[i] = cTors;
1390	–	}
1391	–	}
104		}
105
106	<	void DUFF::fastForward( char* stopText, char* searchOwner ){
1395	<
1396	<	int foundText = 0;
1397	<	char* the_token;
1398	<
1399	<	rewind( frcFile );
1400	<	lineNum = 0;
1401	<
1402	<	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1403	<	lineNum++;
1404	<	if( eof_test == NULL ){
1405	<	sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1406	<	" file is empty.\n",
1407	<	searchOwner );
1408	<	painCave.isFatal = 1;
1409	<	simError();
1410	<	}
1411	<
1412	<
1413	<	while( !foundText ){
1414	<	while( eof_test != NULL && readLine[0] != '#' ){
1415	<	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1416	<	lineNum++;
1417	<	}
1418	<	if( eof_test == NULL ){
1419	<	sprintf( painCave.errMsg,
1420	<	"Error fast forwarding force file for %s at "
1421	<	"line %d: file ended unexpectedly.\n",
1422	<	searchOwner,
1423	<	lineNum );
1424	<	painCave.isFatal = 1;
1425	<	simError();
1426	<	}
1427	<
1428	<	the_token = strtok( readLine, " ,;\t#\n" );
1429	<	foundText = !strcmp( stopText, the_token );
1430	<
1431	<	if( !foundText ){
1432	<	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1433	<	lineNum++;
1434	<
1435	<	if( eof_test == NULL ){
1436	<	sprintf( painCave.errMsg,
1437	<	"Error fast forwarding force file for %s at "
1438	<	"line %d: file ended unexpectedly.\n",
1439	<	searchOwner,
1440	<	lineNum );
1441	<	painCave.isFatal = 1;
1442	<	simError();
1443	<	}
1444	<	}
1445	<	}
1446	<	}
1447	<
1448	<
1449	<	int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1450	<
1451	<	char* the_token;
1452	<
1453	<	the_token = strtok( lineBuffer, " \n\t,;" );
1454	<	if( the_token != NULL ){
1455	<
1456	<	strcpy( info.name, the_token );
1457	<
1458	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1459	<	sprintf( painCave.errMsg,
1460	<	"Error parseing AtomTypes: line %d\n", lineNum );
1461	<	painCave.isFatal = 1;
1462	<	simError();
1463	<	}
1464	<
1465	<	info.isDipole = atoi( the_token );
1466	<
1467	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1468	<	sprintf( painCave.errMsg,
1469	<	"Error parseing AtomTypes: line %d\n", lineNum );
1470	<	painCave.isFatal = 1;
1471	<	simError();
1472	<	}
1473	<
1474	<	info.isSSD = atoi( the_token );
1475	<
1476	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1477	<	sprintf( painCave.errMsg,
1478	<	"Error parseing AtomTypes: line %d\n", lineNum );
1479	<	painCave.isFatal = 1;
1480	<	simError();
1481	<	}
1482	<
1483	<	info.mass = atof( the_token );
1484	<
1485	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1486	<	sprintf( painCave.errMsg,
1487	<	"Error parseing AtomTypes: line %d\n", lineNum );
1488	<	painCave.isFatal = 1;
1489	<	simError();
1490	<	}
1491	<
1492	<	info.epslon = atof( the_token );
1493	<
1494	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1495	<	sprintf( painCave.errMsg,
1496	<	"Error parseing AtomTypes: line %d\n", lineNum );
1497	<	painCave.isFatal = 1;
1498	<	simError();
1499	<	}
1500	<
1501	<	info.sigma = atof( the_token );
1502	<
1503	<	if( info.isDipole ){
1504	<
1505	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1506	<	sprintf( painCave.errMsg,
1507	<	"Error parseing AtomTypes: line %d\n", lineNum );
1508	<	painCave.isFatal = 1;
1509	<	simError();
1510	<	}
1511	<
1512	<	info.dipole = atof( the_token );
1513	<	}
1514	<	else info.dipole = 0.0;
1515	<
1516	<	if( info.isSSD ){
1517	<
1518	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1519	<	sprintf( painCave.errMsg,
1520	<	"Error parseing AtomTypes: line %d\n", lineNum );
1521	<	painCave.isFatal = 1;
1522	<	simError();
1523	<	}
1524	<
1525	<	info.w0 = atof( the_token );
1526	<
1527	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1528	<	sprintf( painCave.errMsg,
1529	<	"Error parseing AtomTypes: line %d\n", lineNum );
1530	<	painCave.isFatal = 1;
1531	<	simError();
1532	<	}
1533	<
1534	<	info.v0 = atof( the_token );
1535	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1536	<	sprintf( painCave.errMsg,
1537	<	"Error parseing AtomTypes: line %d\n", lineNum );
1538	<	painCave.isFatal = 1;
1539	<	simError();
1540	<	}
1541	<
1542	<	info.v0p = atof( the_token );
1543	<
1544	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1545	<	sprintf( painCave.errMsg,
1546	<	"Error parseing AtomTypes: line %d\n", lineNum );
1547	<	painCave.isFatal = 1;
1548	<	simError();
1549	<	}
1550	<
1551	<	info.rl = atof( the_token );
1552	<
1553	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1554	<	sprintf( painCave.errMsg,
1555	<	"Error parseing AtomTypes: line %d\n", lineNum );
1556	<	painCave.isFatal = 1;
1557	<	simError();
1558	<	}
1559	<
1560	<	info.ru = atof( the_token );
1561	<
1562	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1563	<	sprintf( painCave.errMsg,
1564	<	"Error parseing AtomTypes: line %d\n", lineNum );
1565	<	painCave.isFatal = 1;
1566	<	simError();
1567	<	}
1568	<
1569	<	info.rlp = atof( the_token );
1570	<
1571	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1572	<	sprintf( painCave.errMsg,
1573	<	"Error parseing AtomTypes: line %d\n", lineNum );
1574	<	painCave.isFatal = 1;
1575	<	simError();
1576	<	}
1577	<
1578	<	info.rup = atof( the_token );
1579	<	}
1580	<	else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0;
1581	<
1582	<	return 1;
1583	<	}
1584	<	else return 0;
1585	<	}
1586	<
1587	<	int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1588	<
1589	<	char* the_token;
1590	<	char bondType[30];
1591	<
1592	<	the_token = strtok( lineBuffer, " \n\t,;" );
1593	<	if( the_token != NULL ){
1594	<
1595	<	strcpy( info.nameA, the_token );
1596	<
1597	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1598	<	sprintf( painCave.errMsg,
1599	<	"Error parseing BondTypes: line %d\n", lineNum );
1600	<	painCave.isFatal = 1;
1601	<	simError();
1602	<	}
1603	<
1604	<	strcpy( info.nameB, the_token );
1605	<
1606	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1607	<	sprintf( painCave.errMsg,
1608	<	"Error parseing BondTypes: line %d\n", lineNum );
1609	<	painCave.isFatal = 1;
1610	<	simError();
1611	<	}
1612	<
1613	<	strcpy( bondType, the_token );
1614	<
1615	<	if( !strcmp( bondType, "fixed" ) ){
1616	<	info.type = FIXED_BOND;
1617	<
1618	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1619	<	sprintf( painCave.errMsg,
1620	<	"Error parseing BondTypes: line %d\n", lineNum );
1621	<	painCave.isFatal = 1;
1622	<	simError();
1623	<	}
1624	<
1625	<	info.d0 = atof( the_token );
1626	<
1627	<	info.k0=0.0;
1628	<	}
1629	<	else if( !strcmp( bondType, "harmonic" ) ){
1630	<	info.type = HARMONIC_BOND;
1631	<
1632	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1633	<	sprintf( painCave.errMsg,
1634	<	"Error parseing BondTypes: line %d\n", lineNum );
1635	<	painCave.isFatal = 1;
1636	<	simError();
1637	<	}
1638	<
1639	<	info.d0 = atof( the_token );
1640	<
1641	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1642	<	sprintf( painCave.errMsg,
1643	<	"Error parseing BondTypes: line %d\n", lineNum );
1644	<	painCave.isFatal = 1;
1645	<	simError();
1646	<	}
1647	<
1648	<	info.k0 = atof( the_token );
1649	<	}
1650	<
1651	<	else{
1652	<	sprintf( painCave.errMsg,
1653	<	"Unknown DUFF bond type \"%s\" at line %d\n",
1654	<	bondType,
1655	<	lineNum );
1656	<	painCave.isFatal = 1;
1657	<	simError();
1658	<	}
1659	<
1660	<	return 1;
1661	<	}
1662	<	else return 0;
1663	<	}
1664	<
1665	<
1666	<	int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1667	<
1668	<	char* the_token;
1669	<
1670	<	the_token = strtok( lineBuffer, " \n\t,;" );
1671	<	if( the_token != NULL ){
1672	<
1673	<	strcpy( info.nameA, the_token );
1674	<
1675	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1676	<	sprintf( painCave.errMsg,
1677	<	"Error parseing BendTypes: line %d\n", lineNum );
1678	<	painCave.isFatal = 1;
1679	<	simError();
1680	<	}
1681	<
1682	<	strcpy( info.nameB, the_token );
1683	<
1684	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1685	<	sprintf( painCave.errMsg,
1686	<	"Error parseing BendTypes: line %d\n", lineNum );
1687	<	painCave.isFatal = 1;
1688	<	simError();
1689	<	}
1690	<
1691	<	strcpy( info.nameC, the_token );
1692	<
1693	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1694	<	sprintf( painCave.errMsg,
1695	<	"Error parseing BendTypes: line %d\n", lineNum );
1696	<	painCave.isFatal = 1;
1697	<	simError();
1698	<	}
1699	<
1700	<	strcpy( info.type, the_token );
1701	<
1702	<	if( !strcmp( info.type, "quadratic" ) ){
1703	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1704	<	sprintf( painCave.errMsg,
1705	<	"Error parseing BendTypes: line %d\n", lineNum );
1706	<	painCave.isFatal = 1;
1707	<	simError();
1708	<	}
1709	<
1710	<	info.k1 = atof( the_token );
1711	<
1712	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1713	<	sprintf( painCave.errMsg,
1714	<	"Error parseing BendTypes: line %d\n", lineNum );
1715	<	painCave.isFatal = 1;
1716	<	simError();
1717	<	}
1718	<
1719	<	info.k2 = atof( the_token );
1720	<
1721	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1722	<	sprintf( painCave.errMsg,
1723	<	"Error parseing BendTypes: line %d\n", lineNum );
1724	<	painCave.isFatal = 1;
1725	<	simError();
1726	<	}
1727	<
1728	<	info.k3 = atof( the_token );
1729	<
1730	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1731	<	sprintf( painCave.errMsg,
1732	<	"Error parseing BendTypes: line %d\n", lineNum );
1733	<	painCave.isFatal = 1;
1734	<	simError();
1735	<	}
1736	<
1737	<	info.t0 = atof( the_token );
1738	<	}
1739	<
1740	<	else{
1741	<	sprintf( painCave.errMsg,
1742	<	"Unknown DUFF bend type \"%s\" at line %d\n",
1743	<	info.type,
1744	<	lineNum );
1745	<	painCave.isFatal = 1;
1746	<	simError();
1747	<	}
1748	<
1749	<	return 1;
1750	<	}
1751	<	else return 0;
1752	<	}
1753	<
1754	<	int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1755	<
1756	<	char*  the_token;
1757	<
1758	<	the_token = strtok( lineBuffer, " \n\t,;" );
1759	<	if( the_token != NULL ){
1760	<
1761	<	strcpy( info.nameA, the_token );
1762	<
1763	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1764	<	sprintf( painCave.errMsg,
1765	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1766	<	painCave.isFatal = 1;
1767	<	simError();
1768	<	}
1769	<
1770	<	strcpy( info.nameB, the_token );
1771	<
1772	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1773	<	sprintf( painCave.errMsg,
1774	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1775	<	painCave.isFatal = 1;
1776	<	simError();
1777	<	}
1778	<
1779	<	strcpy( info.nameC, the_token );
1780	<
1781	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1782	<	sprintf( painCave.errMsg,
1783	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1784	<	painCave.isFatal = 1;
1785	<	simError();
1786	<	}
1787	<
1788	<	strcpy( info.nameD, the_token );
1789	<
1790	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1791	<	sprintf( painCave.errMsg,
1792	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1793	<	painCave.isFatal = 1;
1794	<	simError();
1795	<	}
1796	<
1797	<	strcpy( info.type, the_token );
1798	<
1799	<	if( !strcmp( info.type, "cubic" ) ){
1800	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1801	<	sprintf( painCave.errMsg,
1802	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1803	<	painCave.isFatal = 1;
1804	<	simError();
1805	<	}
1806	<
1807	<	info.k1 = atof( the_token );
1808	<
1809	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1810	<	sprintf( painCave.errMsg,
1811	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1812	<	painCave.isFatal = 1;
1813	<	simError();
1814	<	}
1815	<
1816	<	info.k2 = atof( the_token );
1817	<
1818	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1819	<	sprintf( painCave.errMsg,
1820	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1821	<	painCave.isFatal = 1;
1822	<	simError();
1823	<	}
1824	<
1825	<	info.k3 = atof( the_token );
1826	<
1827	<	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1828	<	sprintf( painCave.errMsg,
1829	<	"Error parseing TorsionTypes: line %d\n", lineNum );
1830	<	painCave.isFatal = 1;
1831	<	simError();
1832	<	}
1833	<
1834	<	info.k4 = atof( the_token );
1835	<
1836	<	}
1837	<
1838	<	else{
1839	<	sprintf( painCave.errMsg,
1840	<	"Unknown DUFF torsion type \"%s\" at line %d\n",
1841	<	info.type,
1842	<	lineNum );
1843	<	painCave.isFatal = 1;
1844	<	simError();
1845	<	}
1846	<
1847	<	return 1;
1848	<	}
1849	<
1850	<	else return 0;
1851	<	}
106	>	} //end namespace oopse

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines