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

Diff Legend

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