ViewVC Help

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

root/group/branches/new_design/OOPSE-3.0/src/UseTheForce/DUFF.cpp

Comparing:
trunk/OOPSE-3.0/src/UseTheForce/DUFF.cpp (file contents), Revision 1490 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/UseTheForce/DUFF.cpp (file contents), Revision 1759 by tim, Fri Nov 19 18:02:33 2004 UTC

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

Diff Legend

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