ViewVC Help

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

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

Comparing branches/new_design/OOPSE-4/src/UseTheForce/DUFF.cpp (file contents):
Revision 1740 by tim, Wed Nov 3 18:00:36 2004 UTC vs.
Revision 1741 by tim, Tue Nov 16 02:07:14 2004 UTC

#	Line 1 | Line 1
1	<	#include <stdlib.h>
2	<	#include <stdio.h>
3	<	#include <string.h>
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 <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"
26	>	#include "UseTheForce/DUFF.hpp"
27		#include "UseTheForce/DarkSide/lj_interface.h"
13	–	#include "UseTheForce/DarkSide/dipole_interface.h"
14	–	#include "UseTheForce/DarkSide/sticky_interface.h"
28
29	<	#ifdef IS_MPI
17	<	#include "UseTheForce/mpiForceField.h"
18	<	#endif // is_mpi
29	>	namespace oopse {
30
31	+	//definition of createDUFF
32	+	ForceField* createDUFF() {
33	+	return new DUFF();
34	+	}
35
36	<	// define some bond Types
36	>	//register createDUFF to ForceFieldCreator
37	>	ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
38
39	<	#define FIXED_BOND    0
40	<	#define HARMONIC_BOND 1
41	<
42	<
43	<	namespace DUFF_NS {  // restrict the access of the folowing to this file only.
44	<
45	<
30	<	// Declare the structures that will be passed by MPI
31	<
32	<	typedef struct{
33	<	char name[15];
39	>	void DUFF::parse(const std::string& filename) {
40	>	ifstrstream* ffStream;
41	>	ffStream = openForceFiledFile(filename);
42	>	const int bufferSize = 65535;
43	>	char line[bufferSize];
44	>	AtomType* atomType;
45	>	std::string atomTypeName;
46		double mass;
47	<	double epslon;
47	>	double epsilon;
48		double sigma;
49	<	double charge;
50	<	double dipole;
51	<	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	<
49	>	int status;
50	>	int ident = 1; //ident begins from 1 (fortran's index begins from 1)
51	>	int lineNo = 0;
52
53	<	typedef struct{
54	<	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
53	>	while(ffStream.getline(line, bufferSize)){
54	>	++lineNo;
55
56	<	MPI_Datatype mpiAtomStructType;
57	<	MPI_Datatype mpiBondStructType;
58	<	MPI_Datatype mpiBendStructType;
59	<	MPI_Datatype mpiTorsionStructType;
56	>	//a line begins with '#' or '!' is comment which is skipped
57	>	if (line[0] != '#' ||line[0] != '!') {
58	>	StringTokenizer tokenizer(line);
59	>
60	>	if (tokenizer.countTokens() >= 4) {
61	>	atomTypeName = tokenizer.nextToken();
62	>	mass = tokenizer.nextTokenAsFloat();
63	>	epsilon = tokenizer.nextTokenAsFloat();
64	>	sigma = tokenizer.nextTokenAsFloat();
65
66	<	#endif
66	>	atomType = new AtomType();
67	>	atomType->setName(atomTypeName);
68	>	atomType->setMass(mass);
69	>
70	>	//by default, all of the properties in AtomTypeProperties is set to 0
71	>	//In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
72	>	atomType->setLennardJones();
73
74	<	class LinkedAtomType {
75	<	public:
76	<	LinkedAtomType(){
77	<	next = NULL;
78	<	name[0] = '\0';
79	<	}
110	<	~LinkedAtomType(){ if( next != NULL ) delete next; }
74	>	atomType->setIdent(ident);
75	>	atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));
76	>	atomType->addProperty(new DoubleGenericData("Sigma", sigma));
77	>	atomType->complete();
78	>	//notify a new LJtype atom type is created
79	>	newLJtype(&ident, &sigma, &epsilon, &status);
80
81	<	LinkedAtomType* find(const char* key){
82	<	if( !strcmp(name, key) ) return this;
83	<	if( next != NULL ) return next->find(key);
84	<	return NULL;
81	>	//add atom type to AtomTypeContainer
82	>	addAtomType(atomTypeName, atomType);
83	>	++ident;
84	>
85	>	} else {
86	>	sprintf( painCave.errMsg,
87	>	"Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
88	>	"in line %d : %s\n",
89	>	filename.c_str(), lineNo, line);
90	>	painCave.severity = OOPSE_ERROR;
91	>	painCave.isFatal = 1;
92	>	simError();
93	>	}
94	>
95	>	}
96		}
97
98	<	void printMe( void ){
99	<
120	<	std::cerr << "LinkedAtype " << name << ": ident = " << ident << "\n";
121	<	//      if( next != NULL ) next->printMe();
98	>	delete ffStream;
99	>	}
100
123	–	}
101
102	<	void add( atomStruct &info ){
102	>	void DUFF::parseAtomType(){
103	>	}
104
105	<	// check for duplicates
106	<
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	<	}
105	>	void DUFF::parseBondType(){
106	>	}
107
108	<
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(const char* key1, const 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(const char* key1, const char* key2, const 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(const char* key1, const char* key2, const char* key3, const 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
108	>	void DUFF::parseBendType(){
109		}
110
111	+	void DUFF::parseTorsionType(){
112
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	–	at = new DirectionalAtomType();
771	–	else
772	–	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;
113		}
114
115	<
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().c_str() );
1138	<	if( currentAtomType == NULL ){
1139	<	sprintf( painCave.errMsg,
1140	<	"AtomType error, %s not found in force file.\n",
1141	<	the_atoms[i]->getType().c_str() );
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().c_str())){
1164	<	dAtom->setI( waterI );
1165	<	}
1166	<	else if(!strcmp("HEAD",the_atoms[i]->getType().c_str())){
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().c_str() );
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	<
1205	<	Atom** the_atoms;
1206	<	the_atoms = entry_plug->atoms;
1207	<
1208	<
1209	<	// initialize the Bonds
1210	<
1211	<	for( i=0; i<nBonds; i++ ){
1212	<
1213	<	a = the_bonds[i].a;
1214	<	b = the_bonds[i].b;
1215	<
1216	<	currentBondType = headBondType->find(  the_atoms[a]->getType().c_str(),  the_atoms[b]->getType().c_str() );
1217	<	if( currentBondType == NULL ){
1218	<	sprintf( painCave.errMsg,
1219	<	"BondType error, %s - %s not found in force file.\n",
1220	<	the_atoms[a]->getType().c_str(),  the_atoms[b]->getType().c_str() );
1221	<	painCave.isFatal = 1;
1222	<	simError();
1223	<	}
1224	<
1225	<	switch( currentBondType->type ){
1226	<
1227	<	case FIXED_BOND:
1228	<
1229	<	bondArray[i] = new ConstrainedBond( *the_atoms[a],
1230	<	*the_atoms[b],
1231	<	currentBondType->d0 );
1232	<	entry_plug->n_constraints++;
1233	<	break;
1234	<
1235	<	case HARMONIC_BOND:
1236	<
1237	<	bondArray[i] = new HarmonicBond( *the_atoms[a],
1238	<	*the_atoms[b],
1239	<	currentBondType->d0,
1240	<	currentBondType->k0 );
1241	<	break;
1242	<
1243	<	default:
1244	<
1245	<	break;
1246	<	// do nothing
1247	<	}
1248	<	}
1249	<	}
1250	<
1251	<	void DUFF::initializeBends( int nBends, Bend** bendArray,
1252	<	bend_set* the_bends ){
1253	<
1254	<	QuadraticBend* qBend;
1255	<	GhostBend* gBend;
1256	<	Atom** the_atoms;
1257	<	the_atoms = entry_plug->atoms;
1258	<
1259	<	int i, a, b, c;
1260	<	string atomA;
1261	<	string atomB;
1262	<	string atomC;
1263	<
1264	<	// initialize the Bends
1265	<
1266	<	for( i=0; i<nBends; i++ ){
1267	<
1268	<	a = the_bends[i].a;
1269	<	b = the_bends[i].b;
1270	<	c = the_bends[i].c;
1271	<
1272	<	atomA = the_atoms[a]->getType();
1273	<	atomB = the_atoms[b]->getType();
1274	<
1275	<	if( the_bends[i].isGhost ) atomC = "GHOST";
1276	<	else atomC = the_atoms[c]->getType();
1277	<
1278	<	currentBendType = headBendType->find( atomA.c_str(), atomB.c_str(), atomC.c_str());
1279	<	if( currentBendType == NULL ){
1280	<	sprintf( painCave.errMsg, "BendType error, %s - %s - %s not found"
1281	<	" in force file.\n",
1282	<	atomA.c_str(), atomB.c_str(), atomC.c_str() );
1283	<	painCave.isFatal = 1;
1284	<	simError();
1285	<	}
1286	<
1287	<	if( !strcmp( currentBendType->type, "quadratic" ) ){
1288	<
1289	<	if( the_bends[i].isGhost){
1290	<
1291	<	if( the_bends[i].ghost == b ){
1292	<	// do nothing
1293	<	}
1294	<	else if( the_bends[i].ghost == a ){
1295	<	c = a;
1296	<	a = b;
1297	<	b = c;
1298	<	}
1299	<	else{
1300	<	sprintf( painCave.errMsg,
1301	<	"BendType error, %s - %s - %s,\n"
1302	<	"  --> central atom is not "
1303	<	"correctly identified with the "
1304	<	"\"ghostVectorSource = \" tag.\n",
1305	<	atomA.c_str(), atomB.c_str(), atomC.c_str() );
1306	<	painCave.isFatal = 1;
1307	<	simError();
1308	<	}
1309	<
1310	<	gBend = new GhostBend( *the_atoms[a],
1311	<	*the_atoms[b]);
1312	<
1313	<	gBend->setConstants( currentBendType->k1,
1314	<	currentBendType->k2,
1315	<	currentBendType->k3,
1316	<	currentBendType->t0 );
1317	<	bendArray[i] = gBend;
1318	<	}
1319	<	else{
1320	<	qBend = new QuadraticBend( *the_atoms[a],
1321	<	*the_atoms[b],
1322	<	*the_atoms[c] );
1323	<	qBend->setConstants( currentBendType->k1,
1324	<	currentBendType->k2,
1325	<	currentBendType->k3,
1326	<	currentBendType->t0 );
1327	<	bendArray[i] = qBend;
1328	<	}
1329	<	}
1330	<	}
1331	<	}
1332	<
1333	<	void DUFF::initializeTorsions( int nTorsions, Torsion** torsionArray,
1334	<	torsion_set* the_torsions ){
1335	<
1336	<	int i, a, b, c, d;
1337	<
1338	<	CubicTorsion* cTors;
1339	<	Atom** the_atoms;
1340	<	the_atoms = entry_plug->atoms;
1341	<
1342	<	// initialize the Torsions
1343	<
1344	<	for( i=0; i<nTorsions; i++ ){
1345	<
1346	<	a = the_torsions[i].a;
1347	<	b = the_torsions[i].b;
1348	<	c = the_torsions[i].c;
1349	<	d = the_torsions[i].d;
1350	<
1351	<	currentTorsionType = headTorsionType->find( the_atoms[a]->getType().c_str(),
1352	<	the_atoms[b]->getType().c_str(),
1353	<	the_atoms[c]->getType().c_str(),
1354	<	the_atoms[d]->getType().c_str() );
1355	<	if( currentTorsionType == NULL ){
1356	<	sprintf( painCave.errMsg,
1357	<	"TorsionType error, %s - %s - %s - %s not found"
1358	<	" in force file.\n",
1359	<	the_atoms[a]->getType().c_str(),
1360	<	the_atoms[b]->getType().c_str(),
1361	<	the_atoms[c]->getType().c_str(),
1362	<	the_atoms[d]->getType().c_str() );
1363	<	painCave.isFatal = 1;
1364	<	simError();
1365	<	}
1366	<
1367	<	if( !strcmp( currentTorsionType->type, "cubic" ) ){
1368	<
1369	<	cTors = new CubicTorsion( *the_atoms[a], *the_atoms[b],
1370	<	*the_atoms[c], *the_atoms[d] );
1371	<	cTors->setConstants( currentTorsionType->k1, currentTorsionType->k2,
1372	<	currentTorsionType->k3, currentTorsionType->k4 );
1373	<	torsionArray[i] = cTors;
1374	<	}
1375	<	}
1376	<	}
115	>	} //end namespace oopse
116
1378	–	void DUFF::fastForward( char* stopText, char* searchOwner ){
1379	–
1380	–	int foundText = 0;
1381	–	char* the_token;
1382	–
1383	–	rewind( frcFile );
1384	–	lineNum = 0;
1385	–
1386	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1387	–	lineNum++;
1388	–	if( eof_test == NULL ){
1389	–	sprintf( painCave.errMsg, "Error fast forwarding force file for %s: "
1390	–	" file is empty.\n",
1391	–	searchOwner );
1392	–	painCave.isFatal = 1;
1393	–	simError();
1394	–	}
1395	–
1396	–
1397	–	while( !foundText ){
1398	–	while( eof_test != NULL && readLine[0] != '#' ){
1399	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1400	–	lineNum++;
1401	–	}
1402	–	if( eof_test == NULL ){
1403	–	sprintf( painCave.errMsg,
1404	–	"Error fast forwarding force file for %s at "
1405	–	"line %d: file ended unexpectedly.\n",
1406	–	searchOwner,
1407	–	lineNum );
1408	–	painCave.isFatal = 1;
1409	–	simError();
1410	–	}
1411	–
1412	–	the_token = strtok( readLine, " ,;\t#\n" );
1413	–	foundText = !strcmp( stopText, the_token );
1414	–
1415	–	if( !foundText ){
1416	–	eof_test = fgets( readLine, sizeof(readLine), frcFile );
1417	–	lineNum++;
1418	–
1419	–	if( eof_test == NULL ){
1420	–	sprintf( painCave.errMsg,
1421	–	"Error fast forwarding force file for %s at "
1422	–	"line %d: file ended unexpectedly.\n",
1423	–	searchOwner,
1424	–	lineNum );
1425	–	painCave.isFatal = 1;
1426	–	simError();
1427	–	}
1428	–	}
1429	–	}
1430	–	}
1431	–
1432	–
1433	–	int DUFF_NS::parseAtom( char *lineBuffer, int lineNum, atomStruct &info ){
1434	–
1435	–	char* the_token;
1436	–
1437	–	the_token = strtok( lineBuffer, " \n\t,;" );
1438	–	if( the_token != NULL ){
1439	–
1440	–	strcpy( info.name, the_token );
1441	–
1442	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1443	–	sprintf( painCave.errMsg,
1444	–	"Error parseing AtomTypes: line %d\n", lineNum );
1445	–	painCave.isFatal = 1;
1446	–	simError();
1447	–	}
1448	–
1449	–	info.isDipole = atoi( the_token );
1450	–
1451	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1452	–	sprintf( painCave.errMsg,
1453	–	"Error parseing AtomTypes: line %d\n", lineNum );
1454	–	painCave.isFatal = 1;
1455	–	simError();
1456	–	}
1457	–
1458	–	info.isSSD = atoi( the_token );
1459	–
1460	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1461	–	sprintf( painCave.errMsg,
1462	–	"Error parseing AtomTypes: line %d\n", lineNum );
1463	–	painCave.isFatal = 1;
1464	–	simError();
1465	–	}
1466	–
1467	–	info.mass = atof( the_token );
1468	–
1469	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1470	–	sprintf( painCave.errMsg,
1471	–	"Error parseing AtomTypes: line %d\n", lineNum );
1472	–	painCave.isFatal = 1;
1473	–	simError();
1474	–	}
1475	–
1476	–	info.epslon = atof( the_token );
1477	–
1478	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1479	–	sprintf( painCave.errMsg,
1480	–	"Error parseing AtomTypes: line %d\n", lineNum );
1481	–	painCave.isFatal = 1;
1482	–	simError();
1483	–	}
1484	–
1485	–	info.sigma = atof( the_token );
1486	–
1487	–	if( info.isDipole ){
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.dipole = atof( the_token );
1497	–	}
1498	–	else info.dipole = 0.0;
1499	–
1500	–	if( info.isSSD ){
1501	–
1502	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1503	–	sprintf( painCave.errMsg,
1504	–	"Error parseing AtomTypes: line %d\n", lineNum );
1505	–	painCave.isFatal = 1;
1506	–	simError();
1507	–	}
1508	–
1509	–	info.w0 = atof( the_token );
1510	–
1511	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1512	–	sprintf( painCave.errMsg,
1513	–	"Error parseing AtomTypes: line %d\n", lineNum );
1514	–	painCave.isFatal = 1;
1515	–	simError();
1516	–	}
1517	–
1518	–	info.v0 = atof( the_token );
1519	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1520	–	sprintf( painCave.errMsg,
1521	–	"Error parseing AtomTypes: line %d\n", lineNum );
1522	–	painCave.isFatal = 1;
1523	–	simError();
1524	–	}
1525	–
1526	–	info.v0p = atof( the_token );
1527	–
1528	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1529	–	sprintf( painCave.errMsg,
1530	–	"Error parseing AtomTypes: line %d\n", lineNum );
1531	–	painCave.isFatal = 1;
1532	–	simError();
1533	–	}
1534	–
1535	–	info.rl = atof( the_token );
1536	–
1537	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1538	–	sprintf( painCave.errMsg,
1539	–	"Error parseing AtomTypes: line %d\n", lineNum );
1540	–	painCave.isFatal = 1;
1541	–	simError();
1542	–	}
1543	–
1544	–	info.ru = atof( the_token );
1545	–
1546	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1547	–	sprintf( painCave.errMsg,
1548	–	"Error parseing AtomTypes: line %d\n", lineNum );
1549	–	painCave.isFatal = 1;
1550	–	simError();
1551	–	}
1552	–
1553	–	info.rlp = atof( the_token );
1554	–
1555	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1556	–	sprintf( painCave.errMsg,
1557	–	"Error parseing AtomTypes: line %d\n", lineNum );
1558	–	painCave.isFatal = 1;
1559	–	simError();
1560	–	}
1561	–
1562	–	info.rup = atof( the_token );
1563	–	}
1564	–	else info.v0 = info.w0 = info.v0p = info.rl = info.ru = info.rlp = info.rup = 0.0;
1565	–
1566	–	return 1;
1567	–	}
1568	–	else return 0;
1569	–	}
1570	–
1571	–	int DUFF_NS::parseBond( char *lineBuffer, int lineNum, bondStruct &info ){
1572	–
1573	–	char* the_token;
1574	–	char bondType[30];
1575	–
1576	–	the_token = strtok( lineBuffer, " \n\t,;" );
1577	–	if( the_token != NULL ){
1578	–
1579	–	strcpy( info.nameA, the_token );
1580	–
1581	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1582	–	sprintf( painCave.errMsg,
1583	–	"Error parseing BondTypes: line %d\n", lineNum );
1584	–	painCave.isFatal = 1;
1585	–	simError();
1586	–	}
1587	–
1588	–	strcpy( info.nameB, the_token );
1589	–
1590	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1591	–	sprintf( painCave.errMsg,
1592	–	"Error parseing BondTypes: line %d\n", lineNum );
1593	–	painCave.isFatal = 1;
1594	–	simError();
1595	–	}
1596	–
1597	–	strcpy( bondType, the_token );
1598	–
1599	–	if( !strcmp( bondType, "fixed" ) ){
1600	–	info.type = FIXED_BOND;
1601	–
1602	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1603	–	sprintf( painCave.errMsg,
1604	–	"Error parseing BondTypes: line %d\n", lineNum );
1605	–	painCave.isFatal = 1;
1606	–	simError();
1607	–	}
1608	–
1609	–	info.d0 = atof( the_token );
1610	–
1611	–	info.k0=0.0;
1612	–	}
1613	–	else if( !strcmp( bondType, "harmonic" ) ){
1614	–	info.type = HARMONIC_BOND;
1615	–
1616	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1617	–	sprintf( painCave.errMsg,
1618	–	"Error parseing BondTypes: line %d\n", lineNum );
1619	–	painCave.isFatal = 1;
1620	–	simError();
1621	–	}
1622	–
1623	–	info.d0 = atof( the_token );
1624	–
1625	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1626	–	sprintf( painCave.errMsg,
1627	–	"Error parseing BondTypes: line %d\n", lineNum );
1628	–	painCave.isFatal = 1;
1629	–	simError();
1630	–	}
1631	–
1632	–	info.k0 = atof( the_token );
1633	–	}
1634	–
1635	–	else{
1636	–	sprintf( painCave.errMsg,
1637	–	"Unknown DUFF bond type \"%s\" at line %d\n",
1638	–	bondType,
1639	–	lineNum );
1640	–	painCave.isFatal = 1;
1641	–	simError();
1642	–	}
1643	–
1644	–	return 1;
1645	–	}
1646	–	else return 0;
1647	–	}
1648	–
1649	–
1650	–	int DUFF_NS::parseBend( char *lineBuffer, int lineNum, bendStruct &info ){
1651	–
1652	–	char* the_token;
1653	–
1654	–	the_token = strtok( lineBuffer, " \n\t,;" );
1655	–	if( the_token != NULL ){
1656	–
1657	–	strcpy( info.nameA, the_token );
1658	–
1659	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1660	–	sprintf( painCave.errMsg,
1661	–	"Error parseing BendTypes: line %d\n", lineNum );
1662	–	painCave.isFatal = 1;
1663	–	simError();
1664	–	}
1665	–
1666	–	strcpy( info.nameB, the_token );
1667	–
1668	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1669	–	sprintf( painCave.errMsg,
1670	–	"Error parseing BendTypes: line %d\n", lineNum );
1671	–	painCave.isFatal = 1;
1672	–	simError();
1673	–	}
1674	–
1675	–	strcpy( info.nameC, the_token );
1676	–
1677	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1678	–	sprintf( painCave.errMsg,
1679	–	"Error parseing BendTypes: line %d\n", lineNum );
1680	–	painCave.isFatal = 1;
1681	–	simError();
1682	–	}
1683	–
1684	–	strcpy( info.type, the_token );
1685	–
1686	–	if( !strcmp( info.type, "quadratic" ) ){
1687	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1688	–	sprintf( painCave.errMsg,
1689	–	"Error parseing BendTypes: line %d\n", lineNum );
1690	–	painCave.isFatal = 1;
1691	–	simError();
1692	–	}
1693	–
1694	–	info.k1 = atof( the_token );
1695	–
1696	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1697	–	sprintf( painCave.errMsg,
1698	–	"Error parseing BendTypes: line %d\n", lineNum );
1699	–	painCave.isFatal = 1;
1700	–	simError();
1701	–	}
1702	–
1703	–	info.k2 = atof( the_token );
1704	–
1705	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1706	–	sprintf( painCave.errMsg,
1707	–	"Error parseing BendTypes: line %d\n", lineNum );
1708	–	painCave.isFatal = 1;
1709	–	simError();
1710	–	}
1711	–
1712	–	info.k3 = atof( the_token );
1713	–
1714	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1715	–	sprintf( painCave.errMsg,
1716	–	"Error parseing BendTypes: line %d\n", lineNum );
1717	–	painCave.isFatal = 1;
1718	–	simError();
1719	–	}
1720	–
1721	–	info.t0 = atof( the_token );
1722	–	}
1723	–
1724	–	else{
1725	–	sprintf( painCave.errMsg,
1726	–	"Unknown DUFF bend type \"%s\" at line %d\n",
1727	–	info.type,
1728	–	lineNum );
1729	–	painCave.isFatal = 1;
1730	–	simError();
1731	–	}
1732	–
1733	–	return 1;
1734	–	}
1735	–	else return 0;
1736	–	}
1737	–
1738	–	int DUFF_NS::parseTorsion( char *lineBuffer, int lineNum, torsionStruct &info ){
1739	–
1740	–	char*  the_token;
1741	–
1742	–	the_token = strtok( lineBuffer, " \n\t,;" );
1743	–	if( the_token != NULL ){
1744	–
1745	–	strcpy( info.nameA, the_token );
1746	–
1747	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1748	–	sprintf( painCave.errMsg,
1749	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1750	–	painCave.isFatal = 1;
1751	–	simError();
1752	–	}
1753	–
1754	–	strcpy( info.nameB, the_token );
1755	–
1756	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1757	–	sprintf( painCave.errMsg,
1758	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1759	–	painCave.isFatal = 1;
1760	–	simError();
1761	–	}
1762	–
1763	–	strcpy( info.nameC, the_token );
1764	–
1765	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1766	–	sprintf( painCave.errMsg,
1767	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1768	–	painCave.isFatal = 1;
1769	–	simError();
1770	–	}
1771	–
1772	–	strcpy( info.nameD, the_token );
1773	–
1774	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1775	–	sprintf( painCave.errMsg,
1776	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1777	–	painCave.isFatal = 1;
1778	–	simError();
1779	–	}
1780	–
1781	–	strcpy( info.type, the_token );
1782	–
1783	–	if( !strcmp( info.type, "cubic" ) ){
1784	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1785	–	sprintf( painCave.errMsg,
1786	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1787	–	painCave.isFatal = 1;
1788	–	simError();
1789	–	}
1790	–
1791	–	info.k1 = atof( the_token );
1792	–
1793	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1794	–	sprintf( painCave.errMsg,
1795	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1796	–	painCave.isFatal = 1;
1797	–	simError();
1798	–	}
1799	–
1800	–	info.k2 = atof( the_token );
1801	–
1802	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1803	–	sprintf( painCave.errMsg,
1804	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1805	–	painCave.isFatal = 1;
1806	–	simError();
1807	–	}
1808	–
1809	–	info.k3 = atof( the_token );
1810	–
1811	–	if( ( the_token = strtok( NULL, " \n\t,;" ) ) == NULL ){
1812	–	sprintf( painCave.errMsg,
1813	–	"Error parseing TorsionTypes: line %d\n", lineNum );
1814	–	painCave.isFatal = 1;
1815	–	simError();
1816	–	}
1817	–
1818	–	info.k4 = atof( the_token );
1819	–
1820	–	}
1821	–
1822	–	else{
1823	–	sprintf( painCave.errMsg,
1824	–	"Unknown DUFF torsion type \"%s\" at line %d\n",
1825	–	info.type,
1826	–	lineNum );
1827	–	painCave.isFatal = 1;
1828	–	simError();
1829	–	}
1830	–
1831	–	return 1;
1832	–	}
1833	–
1834	–	else return 0;
1835	–	}

Diff Legend

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