ViewVC Help

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

root/OpenMD/branches/devel_omp/src/UseTheForce/DUFF.cpp

Comparing:
trunk/src/UseTheForce/DUFF.cpp (file contents), Revision 143 by chrisfen, Fri Oct 22 22:54:01 2004 UTC vs.
branches/development/src/UseTheForce/DUFF.cpp (file contents), Revision 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC

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

Comparing:
trunk/src/UseTheForce/DUFF.cpp (property svn:keywords), Revision 143 by chrisfen, Fri Oct 22 22:54:01 2004 UTC vs.
branches/development/src/UseTheForce/DUFF.cpp (property svn:keywords), Revision 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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