ViewVC Help

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

root/group/trunk/OOPSE/libmdtools/ZConstraint.cpp

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 660 by tim, Thu Jul 31 19:59:34 2003 UTC vs.
Revision 1203 by gezelter, Thu May 27 18:59:17 2004 UTC

#	Line 1 | Line 1
1		#include "Integrator.hpp"
2		#include "simError.h"
3	+	#include <math.h>
4
5	<	template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
6	<	: T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL)
7	<	{
8	<
5	>	const double INFINITE_TIME = 10e30;
6	>	template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo,
7	>	ForceFields* the_ff): T(theInfo, the_ff),
8	>	fzOut(NULL),
9	>	curZconsTime(0),
10	>	forcePolicy(NULL),
11	>	usingSMD(false),
12	>	hasZConsGap(false){
13		//get properties from SimInfo
14		GenericData* data;
15	<	IndexData* index;
15	>	ZConsParaData* zConsParaData;
16		DoubleData* sampleTime;
17	+	DoubleData* tolerance;
18	+	DoubleData* gap;
19	+	DoubleData* fixtime;
20	+	StringData* policy;
21		StringData* filename;
22	<
23	<
15	<	data = info->getProperty("zconsindex");
16	<	if(!data) {
22	>	IntData* smdFlag;
23	>	double COM[3];
24
25	<	sprintf( painCave.errMsg,
26	<	"ZConstraint error: If you use an ZConstraint\n"
27	<	" , you must set index of z-constraint molecules.\n");
28	<	painCave.isFatal = 1;
29	<	simError();
25	>	//by default, the direction of constraint is z
26	>	// 0 --> x
27	>	// 1 --> y
28	>	// 2 --> z
29	>	whichDirection = 2;
30	>
31	>	//estimate the force constant of harmonical potential
32	>	double Kb = 1.986E-3 ; //in kcal/K
33	>
34	>	double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /
35	>	2;
36	>	zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox);
37	>
38	>	//creat force Subtraction policy
39	>	data = info->getProperty(ZCONSFORCEPOLICY_ID);
40	>	if (!data){
41	>	sprintf(painCave.errMsg,
42	>	"ZConstraint Warning: User does not set force Subtraction policy, "
43	>	"PolicyByMass is used\n");
44	>	painCave.isFatal = 0;
45	>	simError();
46	>
47	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
48		}
49		else{
50	<	index = dynamic_cast<IndexData*>(data);
26	<
27	<	if(!index){
50	>	policy = dynamic_cast<StringData*>(data);
51
52	<	sprintf( painCave.errMsg,
53	<	"ZConstraint error: Can not get property from SimInfo\n");
54	<	painCave.isFatal = 1;
55	<	simError();
56	<
52	>	if (!policy){
53	>	sprintf(painCave.errMsg,
54	>	"ZConstraint Error: Convertion from GenericData to StringData failure, "
55	>	"PolicyByMass is used\n");
56	>	painCave.isFatal = 0;
57	>	simError();
58	>
59	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
60		}
61		else{
62	<
63	<	indexOfAllZConsMols = index->getIndexData();
64	<
65	<	//the maximum value of index is the last one(we sorted the index data in SimSetup.cpp)
66	<	int maxIndex;
67	<	int totalNumMol;
68	<
69	<	maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1];
70	<
71	<	#ifndef IS_MPI
72	<	totalNumMol = nMols;
47	<	#else
48	<	totalNumMol = mpiSim->getTotNmol();
49	<	#endif
50	<
51	<	if(maxIndex > totalNumMol - 1){
52	<	sprintf( painCave.errMsg,
53	<	"ZConstraint error: index is out of range\n");
54	<	painCave.isFatal = 1;
55	<	simError();
56	<
62	>	if (policy->getData() == "BYNUMBER")
63	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this);
64	>	else if (policy->getData() == "BYMASS")
65	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
66	>	else{
67	>	sprintf(painCave.errMsg,
68	>	"ZConstraint Warning: unknown force Subtraction policy, "
69	>	"PolicyByMass is used\n");
70	>	painCave.isFatal = 0;
71	>	simError();
72	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
73		}
58	–
74		}
60	–
75		}
76	<
77	<	//retrive sample time of z-contraint
78	<	data = info->getProperty("zconstime");
79	<
80	<	if(!data) {
81	<
82	<	sprintf( painCave.errMsg,
83	<	"ZConstraint error: If you use an ZConstraint\n"
84	<	" , you must set sample time.\n");
76	>
77	>
78	>	//retrieve sample time of z-contraint
79	>	data = info->getProperty(ZCONSTIME_ID);
80	>
81	>	if (!data){
82	>	sprintf(painCave.errMsg,
83	>	"ZConstraint error: If you use an ZConstraint\n"
84	>	" , you must set sample time.\n");
85		painCave.isFatal = 1;
86	<	simError();
86	>	simError();
87		}
88		else{
75	–
89		sampleTime = dynamic_cast<DoubleData*>(data);
77	–
78	–	if(!sampleTime){
90
91	<	sprintf( painCave.errMsg,
92	<	"ZConstraint error: Can not get property from SimInfo\n");
91	>	if (!sampleTime){
92	>	sprintf(painCave.errMsg,
93	>	"ZConstraint error: Can not get property from SimInfo\n");
94		painCave.isFatal = 1;
95	<	simError();
84	<
95	>	simError();
96		}
97		else{
98		this->zconsTime = sampleTime->getData();
99		}
89	–
100		}
91	–
92	–
93	–	//retrive output filename of z force
94	–	data = info->getProperty("zconsfilename");
95	–	if(!data) {
101
102	<
103	<	sprintf( painCave.errMsg,
104	<	"ZConstraint error: If you use an ZConstraint\n"
105	<	" , you must set output filename of z-force.\n");
102	>	//retrieve output filename of z force
103	>	data = info->getProperty(ZCONSFILENAME_ID);
104	>	if (!data){
105	>	sprintf(painCave.errMsg,
106	>	"ZConstraint error: If you use an ZConstraint\n"
107	>	" , you must set output filename of z-force.\n");
108		painCave.isFatal = 1;
109	<	simError();
103	<
109	>	simError();
110		}
111		else{
112	+	filename = dynamic_cast<StringData*>(data);
113
114	<	filename = dynamic_cast<StringData*>(data);
115	<
116	<	if(!filename){
110	<
111	<	sprintf( painCave.errMsg,
112	<	"ZConstraint error: Can not get property from SimInfo\n");
114	>	if (!filename){
115	>	sprintf(painCave.errMsg,
116	>	"ZConstraint error: Can not get property from SimInfo\n");
117		painCave.isFatal = 1;
118	<	simError();
115	<
118	>	simError();
119		}
120		else{
121		this->zconsOutput = filename->getData();
122		}
120	–
121	–
123		}
124
125	<
126	<	//calculate reference z coordinate for z-constraint molecules
126	<	double totalMass_local;
127	<	double totalMass;
128	<	double totalMZ_local;
129	<	double totalMZ;
130	<	double massOfUncons_local;
131	<	double massOfCurMol;
132	<	double COM[3];
133	<
134	<	totalMass_local = 0;
135	<	totalMass = 0;
136	<	totalMZ_local = 0;
137	<	totalMZ = 0;
138	<	massOfUncons_local = 0;
139	<
140	<
141	<	for(int i = 0; i < nMols; i++){
142	<	massOfCurMol = molecules[i].getTotalMass();
143	<	molecules[i].getCOM(COM);
144	<
145	<	totalMass_local += massOfCurMol;
146	<	totalMZ_local += massOfCurMol * COM[2];
147	<
148	<	if(isZConstraintMol(&molecules[i]) == -1){
149	<
150	<	massOfUncons_local += massOfCurMol;
151	<	}
152	<
153	<	}
154	<
155	<
156	<	#ifdef IS_MPI
157	<	MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
158	<	MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
159	<	MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
160	<	#else
161	<	totalMass = totalMass_local;
162	<	totalMZ = totalMZ_local;
163	<	totalMassOfUncons = massOfUncons_local;
164	<	#endif
125	>	//retrieve tolerance for z-constraint molecuels
126	>	data = info->getProperty(ZCONSTOL_ID);
127
128	<	double zsys;
129	<	zsys = totalMZ / totalMass;
128	>	if (!data){
129	>	sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n");
130	>	painCave.isFatal = 1;
131	>	simError();
132	>	}
133	>	else{
134	>	tolerance = dynamic_cast<DoubleData*>(data);
135
136	<	#ifndef IS_MPI
137	<	for(int i = 0; i < nMols; i++){
138	<
139	<	if(isZConstraintMol(&molecules[i]) > -1 ){
140	<	molecules[i].getCOM(COM);
174	<	allRefZ.push_back(COM[2] - zsys);
136	>	if (!tolerance){
137	>	sprintf(painCave.errMsg,
138	>	"ZConstraint error: Can not get property from SimInfo\n");
139	>	painCave.isFatal = 1;
140	>	simError();
141		}
142	<
142	>	else{
143	>	this->zconsTol = tolerance->getData();
144	>	}
145		}
146	<	#else
147	<
148	<	int whichNode;
149	<	enum CommType { RequestMolZPos, EndOfRequest} status;
150	<	//int status;
151	<	double zpos;
152	<	int localIndex;
153	<	MPI_Status ierr;
154	<	int tag = 0;
155	<
156	<	if(worldRank == 0){
157	<
158	<	int globalIndexOfCurMol;
159	<	int *MolToProcMap;
160	<	MolToProcMap = mpiSim->getMolToProcMap();
161	<
162	<	for(int i = 0; i < indexOfAllZConsMols.size(); i++){
195	<
196	<	whichNode = MolToProcMap[indexOfAllZConsMols[i]];
197	<	globalIndexOfCurMol = indexOfAllZConsMols[i];
198	<
199	<	if(whichNode == 0){
200	<
201	<	for(int j = 0; j < nMols; j++)
202	<	if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){
203	<	localIndex = j;
204	<	break;
205	<	}
206	<
207	<	molecules[localIndex].getCOM(COM);
208	<	allRefZ.push_back(COM[2] - zsys);
209	<
210	<	}
211	<	else{
212	<	status = RequestMolZPos;
213	<	MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
214	<	MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
215	<	MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr);
216	<
217	<	allRefZ.push_back(zpos - zsys);
218	<
219	<	}
220	<
221	<	} //End of Request Loop
222	<
223	<	//Send ending request message to slave nodes
224	<	status = EndOfRequest;
225	<	for(int i =1; i < mpiSim->getNumberProcessors(); i++)
226	<	MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
227	<
146	>
147	>	//quick hack here
148	>	data = info->getProperty(ZCONSGAP_ID);
149	>
150	>	if (data){
151	>	gap = dynamic_cast<DoubleData*>(data);
152	>
153	>	if (!gap){
154	>	sprintf(painCave.errMsg,
155	>	"ZConstraint error: Can not get property from SimInfo\n");
156	>	painCave.isFatal = 1;
157	>	simError();
158	>	}
159	>	else{
160	>	this->hasZConsGap = true;
161	>	this->zconsGap = gap->getData();
162	>	}
163		}
229	–	else{
230	–
231	–	int whichMol;
232	–	bool done = false;
164
165	<	while (!done){
166	<
167	<	MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
168	<
169	<	switch (status){
170	<
171	<	case RequestMolZPos :
172	<
173	<	MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr);
174	<
175	<	for(int i = 0; i < nMols; i++)
245	<	if(molecules[i].getGlobalIndex() == whichMol){
246	<	localIndex = i;
247	<	break;
248	<	}
249	<
250	<	molecules[localIndex].getCOM(COM);
251	<	zpos = COM[2];
252	<	MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD);
253	<
254	<	break;
255	<
256	<	case EndOfRequest :
257	<
258	<	done = true;
259	<	break;
260	<	}
261	<
165	>
166	>
167	>	data = info->getProperty(ZCONSFIXTIME_ID);
168	>
169	>	if (data){
170	>	fixtime = dynamic_cast<DoubleData*>(data);
171	>	if (!fixtime){
172	>	sprintf(painCave.errMsg,
173	>	"ZConstraint error: Can not get zconsFixTime from SimInfo\n");
174	>	painCave.isFatal = 1;
175	>	simError();
176		}
177	<
177	>	else{
178	>	this->zconsFixTime = fixtime->getData();
179	>	}
180		}
181	+	else if(hasZConsGap){
182	+	sprintf(painCave.errMsg,
183	+	"ZConstraint error: must set fixtime if already set zconsGap\n");
184	+	painCave.isFatal = 1;
185	+	simError();
186	+	}
187
266	–	//Brocast the allRefZ to slave nodes;
267	–	double* allRefZBuf;
268	–	int nZConsMols;
269	–	nZConsMols = indexOfAllZConsMols.size();
270	–
271	–	allRefZBuf = new double[nZConsMols];
272	–
273	–	if(worldRank == 0){
188
189	<	for(int i = 0; i < nZConsMols; i++)
190	<	allRefZBuf[i] = allRefZ[i];
191	<	}
192	<
193	<	MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
194	<
195	<	if(worldRank != 0){
196	<
197	<	for(int i = 0; i < nZConsMols; i++)
198	<	allRefZ.push_back(allRefZBuf[i]);
189	>
190	>	data = info->getProperty(ZCONSUSINGSMD_ID);
191	>
192	>	if (data){
193	>	smdFlag = dynamic_cast<IntData*>(data);
194	>
195	>	if (!smdFlag){
196	>	sprintf(painCave.errMsg,
197	>	"ZConstraint error: Can not get property from SimInfo\n");
198	>	painCave.isFatal = 1;
199	>	simError();
200	>	}
201	>	else{
202	>	this->usingSMD= smdFlag->getData() ? true : false;
203	>	}
204	>
205		}
206	<
207	<	delete[] allRefZBuf;
206	>
207	>
208	>
209	>	//retrieve index of z-constraint molecules
210	>	data = info->getProperty(ZCONSPARADATA_ID);
211	>	if (!data){
212	>	sprintf(painCave.errMsg,
213	>	"ZConstraint error: If you use an ZConstraint\n"
214	>	" , you must set index of z-constraint molecules.\n");
215	>	painCave.isFatal = 1;
216	>	simError();
217	>	}
218	>	else{
219	>	zConsParaData = dynamic_cast<ZConsParaData*>(data);
220	>
221	>	if (!zConsParaData){
222	>	sprintf(painCave.errMsg,
223	>	"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
224	>	painCave.isFatal = 1;
225	>	simError();
226	>	}
227	>	else{
228	>	parameters = zConsParaData->getData();
229	>
230	>	//check the range of zconsIndex
231	>	//and the minimum value of index is the first one (we already sorted the data)
232	>	//the maximum value of index is the last one
233	>
234	>	int maxIndex;
235	>	int minIndex;
236	>	int totalNumMol;
237	>
238	>	minIndex = (*parameters)[0].zconsIndex;
239	>	if (minIndex < 0){
240	>	sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
241	>	painCave.isFatal = 1;
242	>	simError();
243	>	}
244	>
245	>	maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
246	>
247	>	#ifndef IS_MPI
248	>	totalNumMol = nMols;
249	>	#else
250	>	totalNumMol = mpiSim->getNMolGlobal();
251	>	#endif
252	>
253	>	if (maxIndex > totalNumMol - 1){
254	>	sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
255	>	painCave.isFatal = 1;
256	>	simError();
257	>	}
258	>
259	>	//if user does not specify the zpos for the zconstraint molecule
260	>	//its initial z coordinate  will be used as default
261	>	for (int i = 0; i < (int) (parameters->size()); i++){
262	>	if (!(*parameters)[i].havingZPos){
263	>	#ifndef IS_MPI
264	>	for (int j = 0; j < nMols; j++){
265	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
266	>	molecules[j].getCOM(COM);
267	>	break;
268	>	}
269	>	}
270	>	#else
271	>	//query which processor current zconstraint molecule belongs to
272	>	int* MolToProcMap;
273	>	int whichNode;
274	>
275	>	MolToProcMap = mpiSim->getMolToProcMap();
276	>	whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
277	>
278	>	//broadcast the zpos of current z-contraint molecule
279	>	//the node which contain this
280	>
281	>	if (worldRank == whichNode){
282	>	for (int j = 0; j < nMols; j++)
283	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
284	>	molecules[j].getCOM(COM);
285	>	break;
286	>	}
287	>	}
288	>
289	>	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode,
290	>	MPI_COMM_WORLD);
291		#endif
292
293	<
293	>	(*parameters)[i].zPos = COM[whichDirection];
294	>
295	>	sprintf(painCave.errMsg,
296	>	"ZConstraint warning: Does not specify zpos for z-constraint molecule "
297	>	"initial z coornidate will be used \n");
298	>	painCave.isFatal = 0;
299	>	simError();
300	>	}
301	>	}
302	>	}//end if (!zConsParaData)
303	>
304	>	}//end  if (!data)
305	>
306	>	//
307		#ifdef IS_MPI
308		update();
309		#else
310		int searchResult;
295	–
296	–	refZ = allRefZ;
311
312	<	for(int i = 0; i < nMols; i++){
299	<
312	>	for (int i = 0; i < nMols; i++){
313		searchResult = isZConstraintMol(&molecules[i]);
314	<
315	<	if(searchResult > -1){
303	<
314	>
315	>	if (searchResult > -1){
316		zconsMols.push_back(&molecules[i]);
317		massOfZConsMols.push_back(molecules[i].getTotalMass());
318	<
319	<	molecules[i].getCOM(COM);
318	>
319	>	zPos.push_back((*parameters)[searchResult].zPos);
320	>	kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
321	>
322	>	if(usingSMD)
323	>	cantVel.push_back((*parameters)[searchResult].cantVel);
324	>
325		}
326	<	else
310	<	{
311	<
326	>	else{
327		unconsMols.push_back(&molecules[i]);
328		massOfUnconsMols.push_back(molecules[i].getTotalMass());
314	–
329		}
330		}
331
332	<	fz = new double[zconsMols.size()];
333	<	indexOfZConsMols = new int [zconsMols.size()];
332	>	fz.resize(zconsMols.size());
333	>	curZPos.resize(zconsMols.size());
334	>	indexOfZConsMols.resize(zconsMols.size());
335
336	<	if(!fz || !indexOfZConsMols){
337	<	sprintf( painCave.errMsg,
338	<	"Memory allocation failure in class Zconstraint\n");
339	<	painCave.isFatal = 1;
340	<	simError();
336	>	//determine the states of z-constraint molecules
337	>	for (size_t i = 0; i < zconsMols.size(); i++){
338	>	indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
339	>
340	>	zconsMols[i]->getCOM(COM);
341	>
342	>	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
343	>	states.push_back(zcsFixed);
344	>
345	>	if (hasZConsGap)
346	>	endFixTime.push_back(info->getTime() + zconsFixTime);
347	>	}
348	>	else{
349	>	states.push_back(zcsMoving);
350	>
351	>	if (hasZConsGap)
352	>	endFixTime.push_back(INFINITE_TIME);
353	>	}
354	>
355	>	if(usingSMD)
356	>	cantPos.push_back(COM[whichDirection]);
357		}
358
359	<	for(int i = 0; i < zconsMols.size(); i++)
360	<	indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
330	<
359	>	if(usingSMD)
360	>	prevCantPos = cantPos;
361		#endif
362	+
363
364	<	fzOut = new ZConsWriter(zconsOutput.c_str());
365	<
366	<	if(!fzOut){
367	<	sprintf( painCave.errMsg,
368	<	"Memory allocation failure in class Zconstraint\n");
369	<	painCave.isFatal = 1;
370	<	simError();
371	<	}
372	<
373	<	fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
364	>	//get total masss of unconstraint molecules
365	>	double totalMassOfUncons_local;
366	>	totalMassOfUncons_local = 0;
367	>
368	>	for (size_t i = 0; i < unconsMols.size(); i++)
369	>	totalMassOfUncons_local += unconsMols[i]->getTotalMass();
370	>
371	>	#ifndef IS_MPI
372	>	totalMassOfUncons = totalMassOfUncons_local;
373	>	#else
374	>	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,
375	>	MPI_SUM, MPI_COMM_WORLD);
376	>	#endif
377	>
378	>	//get total number of unconstrained atoms
379	>	int nUnconsAtoms_local;
380	>	nUnconsAtoms_local = 0;
381	>	for (int i = 0; i < (int) (unconsMols.size()); i++)
382	>	nUnconsAtoms_local += unconsMols[i]->getNAtoms();
383	>
384	>	#ifndef IS_MPI
385	>	totNumOfUnconsAtoms = nUnconsAtoms_local;
386	>	#else
387	>	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM,
388	>	MPI_COMM_WORLD);
389	>	#endif
390	>
391	>	forcePolicy->update();
392		}
393
394	<	template<typename T> ZConstraint<T>::~ZConstraint()
395	<	{
396	<	if(fz)
348	<	delete[] fz;
349	<
350	<	if(indexOfZConsMols)
351	<	delete[] indexOfZConsMols;
352	<
353	<	if(fzOut)
394	>	template<typename T> ZConstraint<T>::~ZConstraint(){
395	>
396	>	if (fzOut){
397		delete fzOut;
398	+	}
399	+
400	+	if (forcePolicy){
401	+	delete forcePolicy;
402	+	}
403		}
404
405	+
406	+	/**
407	+	*
408	+	*/
409	+
410		#ifdef IS_MPI
411	<	template<typename T> void ZConstraint<T>::update()
359	<	{
411	>	template<typename T> void ZConstraint<T>::update(){
412		double COM[3];
413		int index;
414	<
414	>
415		zconsMols.clear();
416		massOfZConsMols.clear();
417	<	refZ.clear();
418	<
417	>	zPos.clear();
418	>	kz.clear();
419	>	cantPos.clear();
420	>	cantVel.clear();
421	>
422		unconsMols.clear();
423		massOfUnconsMols.clear();
369	–
424
425	+
426		//creat zconsMol and unconsMol lists
427	<	for(int i = 0; i < nMols; i++){
373	<
427	>	for (int i = 0; i < nMols; i++){
428		index = isZConstraintMol(&molecules[i]);
429	<
430	<	if(index > -1){
377	<
429	>
430	>	if (index > -1){
431		zconsMols.push_back(&molecules[i]);
432	<	massOfZConsMols.push_back(molecules[i].getTotalMass());
433	<
434	<	molecules[i].getCOM(COM);
435	<	refZ.push_back(allRefZ[index]);
432	>	zPos.push_back((*parameters)[index].zPos);
433	>	kz.push_back((*parameters)[index].kRatio * zForceConst);
434	>	massOfZConsMols.push_back(molecules[i].getTotalMass());
435	>
436	>	if(usingSMD)
437	>	cantVel.push_back((*parameters)[index].cantVel);
438	>
439		}
440	<	else
385	<	{
386	<
440	>	else{
441		unconsMols.push_back(&molecules[i]);
442		massOfUnconsMols.push_back(molecules[i].getTotalMass());
389	–
443		}
444		}
445	+
446	+	fz.resize(zconsMols.size());
447	+	curZPos.resize(zconsMols.size());
448	+	indexOfZConsMols.resize(zconsMols.size());
449	+
450	+	for (size_t i = 0; i < zconsMols.size(); i++){
451	+	indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
452	+	}
453
454	<	//The reason to declare fz and indexOfZconsMols as pointer to array is
455	<	// that we want to make the MPI communication simple
456	<	if(fz)
457	<	delete[] fz;
454	>	//determine the states of z-constraint molecules
455	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
456	>
457	>	zconsMols[i]->getCOM(COM);
458
459	<	if(indexOfZConsMols)
460	<	delete[] indexOfZConsMols;
461	<
462	<	if (zconsMols.size() > 0){
463	<	fz = new double[zconsMols.size()];
403	<	indexOfZConsMols =  new int[zconsMols.size()];
404	<
405	<	if(!fz || !indexOfZConsMols){
406	<	sprintf( painCave.errMsg,
407	<	"Memory allocation failure in class Zconstraint\n");
408	<	painCave.isFatal = 1;
409	<	simError();
459	>	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
460	>	states.push_back(zcsFixed);
461	>
462	>	if (hasZConsGap)
463	>	endFixTime.push_back(info->getTime() + zconsFixTime);
464		}
465	<
466	<	for(int i = 0; i < zconsMols.size(); i++){
467	<	indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
465	>	else{
466	>	states.push_back(zcsMoving);
467	>
468	>	if (hasZConsGap)
469	>	endFixTime.push_back(INFINITE_TIME);
470		}
471
472	<	}
473	<	else{
474	<	fz = NULL;
475	<	indexOfZConsMols = NULL;
476	<	}
477	<
472	>	if(usingSMD)
473	>	cantPos.push_back(COM[whichDirection]);
474	>	}
475	>
476	>	if(usingSMD)
477	>	prevCantPos = cantPos;
478	>
479	>	//
480	>	forcePolicy->update();
481		}
482
483		#endif
484
485	<	/**  Function Name: isZConstraintMol
486	<	**  Parameter
487	<	**    Molecule* mol
488	<	**  Return value:
489	<	**    -1, if the molecule is not z-constraint molecule,
490	<	**    other non-negative values, its index in indexOfAllZConsMols vector
485	>	/**
486	>	*  Function Name: isZConstraintMol
487	>	*  Parameter
488	>	*    Molecule* mol
489	>	*  Return value:
490	>	*    -1, if the molecule is not z-constraint molecule,
491	>	*    other non-negative values, its index in indexOfAllZConsMols vector
492		*/
493
494	<	template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
435	<	{
494	>	template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
495		int index;
496		int low;
497		int high;
498		int mid;
499
500		index = mol->getGlobalIndex();
501	<
501	>
502		low = 0;
503	<	high = indexOfAllZConsMols.size() - 1;
504	<
503	>	high = parameters->size() - 1;
504	>
505		//Binary Search (we have sorted the array)
506	<	while(low <= high){
507	<	mid = (low + high) /2;
508	<	if (indexOfAllZConsMols[mid] == index)
506	>	while (low <= high){
507	>	mid = (low + high) / 2;
508	>	if ((*parameters)[mid].zconsIndex == index)
509		return mid;
510	<	else if (indexOfAllZConsMols[mid] > index )
511	<	high = mid -1;
512	<	else
513	<	low = mid + 1;
510	>	else if ((*parameters)[mid].zconsIndex > index)
511	>	high = mid - 1;
512	>	else
513	>	low = mid + 1;
514		}
515	<
515	>
516		return -1;
517		}
518
519	<	/** Function Name: integrateStep
520	<	** Parameter:
521	<	**   int calcPot;
522	<	**   int calcStress;
523	<	** Description:
524	<	**  Advance One Step.
525	<	** Memo:
526	<	**   The best way to implement z-constraint is to override integrateStep
527	<	**   Overriding constrainB is not a good choice, since in integrateStep,
528	<	**   constrainB is invoked by below line,
529	<	**                  if(nConstrained) constrainB();
530	<	**   For instance, we would like to apply z-constraint without bond contrain,
531	<	**   In that case, if we override constrainB, Z-constrain method will never be executed;
532	<	*/
533	<	template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
534	<	{
535	<	T::integrateStep( calcPot, calcStress );
477	<	resetZ();
478	<
479	<	double currZConsTime = 0;
480	<
481	<	//write out forces of z constraint
482	<	if( info->getTime() >= currZConsTime){
483	<	fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
484	<	}
519	>	template<typename T> void ZConstraint<T>::integrate(){
520	>	// creat zconsWriter
521	>	fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
522	>
523	>	if (!fzOut){
524	>	sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
525	>	painCave.isFatal = 1;
526	>	simError();
527	>	}
528	>
529	>	//zero out the velocities of center of mass of unconstrained molecules
530	>	//and the velocities of center of mass of every single z-constrained molecueles
531	>	zeroOutVel();
532	>
533	>	curZconsTime = zconsTime + info->getTime();
534	>
535	>	T::integrate();
536		}
537
538	<	/** Function Name: resetZ
539	<	** Description:
540	<	**  Reset the z coordinates
538	>
539	>	/**
540	>	*
541	>	*
542	>	*
543	>	*
544	>	*/
545	>	template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
546	>	double zsys;
547	>	double COM[3];
548	>	double force[3];
549	>	double zSysCOMVel;
550	>
551	>	T::calcForce(calcPot, calcStress);
552	>
553	>
554	>	if (hasZConsGap){
555	>	updateZPos();
556	>	}
557	>
558	>	if (checkZConsState()){
559	>	zeroOutVel();
560	>	forcePolicy->update();
561	>	}
562	>
563	>	zsys = calcZSys();
564	>	zSysCOMVel = calcSysCOMVel();
565	>	#ifdef IS_MPI
566	>	if (worldRank == 0){
567	>	#endif
568	>	//cout << "---------------------------------------------------------------------" <<endl;
569	>	//cout << "current time: " << info->getTime() << endl;
570	>	//cout << "center of mass at z: " << zsys << endl;
571	>	//cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
572	>
573	>	#ifdef IS_MPI
574	>	}
575	>	#endif
576	>
577	>	//do zconstraint force;
578	>	if (haveFixedZMols()){
579	>	this->doZconstraintForce();
580	>	}
581	>
582	>	//use external force to move the molecules to the specified positions
583	>	if (haveMovingZMols()){
584	>	if (usingSMD)
585	>	this->doHarmonic(cantPos);
586	>	else
587	>	this->doHarmonic(zPos);
588	>	}
589	>
590	>	//write out forces and current positions of z-constraint molecules
591	>	if (info->getTime() >= curZconsTime){
592	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
593	>	zconsMols[i]->getCOM(COM);
594	>	curZPos[i] = COM[whichDirection];
595	>
596	>	//if the z-constraint molecule is still moving, just record its force
597	>	if (states[i] == zcsMoving){
598	>	fz[i] = 0;
599	>	Atom** movingZAtoms;
600	>	movingZAtoms = zconsMols[i]->getMyAtoms();
601	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
602	>	movingZAtoms[j]->getFrc(force);
603	>	fz[i] += force[whichDirection];
604	>	}
605	>	}
606	>	}
607	>	fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
608	>	&curZPos[0], &zPos[0]);
609	>	curZconsTime += zconsTime;
610	>	}
611	>
612	>	zSysCOMVel = calcSysCOMVel();
613	>	#ifdef IS_MPI
614	>	if (worldRank == 0){
615	>	#endif
616	>	//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
617	>	#ifdef IS_MPI
618	>	}
619	>	#endif
620	>	}
621	>
622	>
623	>	/**
624	>	*
625		*/
626
627	<	template<typename T> void ZConstraint<T>::resetZ()
628	<	{
629	<	double deltaZ;
630	<	double mzOfZCons;   //total sum of m*z of z-constrain molecules
631	<	double mzOfUncons; //total sum of m*z of unconstrain molecuels;
632	<	double totalMZOfZCons;
633	<	double totalMZOfUncons;
627	>	template<typename T> double ZConstraint<T>::calcZSys(){
628	>	//calculate reference z coordinate for z-constraint molecules
629	>	double totalMass_local;
630	>	double totalMass;
631	>	double totalMZ_local;
632	>	double totalMZ;
633	>	double massOfCurMol;
634		double COM[3];
635	+
636	+	totalMass_local = 0;
637	+	totalMZ_local = 0;
638	+
639	+	for (int i = 0; i < nMols; i++){
640	+	massOfCurMol = molecules[i].getTotalMass();
641	+	molecules[i].getCOM(COM);
642	+
643	+	totalMass_local += massOfCurMol;
644	+	totalMZ_local += massOfCurMol * COM[whichDirection];
645	+	}
646	+
647	+
648	+	#ifdef IS_MPI
649	+	MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
650	+	MPI_COMM_WORLD);
651	+	MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
652	+	#else
653	+	totalMass = totalMass_local;
654	+	totalMZ = totalMZ_local;
655	+	#endif
656	+
657		double zsys;
658	<	Atom** zconsAtoms;
658	>	zsys = totalMZ / totalMass;
659
660	<	mzOfZCons = 0;
661	<	mzOfUncons  = 0;
662	<
663	<	for(int i = 0; i < zconsMols.size(); i++){
664	<	mzOfZCons += massOfZConsMols[i] * refZ[i];
660	>	return zsys;
661	>	}
662	>
663	>	/**
664	>	*
665	>	*/
666	>	template<typename T> void ZConstraint<T>::thermalize(void){
667	>	T::thermalize();
668	>	zeroOutVel();
669	>	}
670	>
671	>	/**
672	>	*
673	>	*/
674	>
675	>	template<typename T> void ZConstraint<T>::zeroOutVel(){
676	>	Atom** fixedZAtoms;
677	>	double COMvel[3];
678	>	double vel[3];
679	>	double zSysCOMVel;
680	>
681	>	//zero out the velocities of center of mass of fixed z-constrained molecules
682	>
683	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
684	>	if (states[i] == zcsFixed){
685	>	zconsMols[i]->getCOMvel(COMvel);
686	>	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
687	>
688	>	fixedZAtoms = zconsMols[i]->getMyAtoms();
689	>
690	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
691	>	fixedZAtoms[j]->getVel(vel);
692	>	vel[whichDirection] -= COMvel[whichDirection];
693	>	fixedZAtoms[j]->setVel(vel);
694	>	}
695	>
696	>	zconsMols[i]->getCOMvel(COMvel);
697	>	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
698	>	}
699		}
700
701	+	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
702	+
703	+	zSysCOMVel = calcSysCOMVel();
704		#ifdef IS_MPI
705	<	MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
705	>	if (worldRank == 0){
706	>	#endif
707	>	//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
708	>	#ifdef IS_MPI
709	>	}
710	>	#endif
711	>
712	>	// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
713	>	double MVzOfMovingMols_local;
714	>	double MVzOfMovingMols;
715	>	double totalMassOfMovingZMols_local;
716	>	double totalMassOfMovingZMols;
717	>
718	>	MVzOfMovingMols_local = 0;
719	>	totalMassOfMovingZMols_local = 0;
720	>
721	>	for (int i = 0; i < (int) (unconsMols.size()); i++){
722	>	unconsMols[i]->getCOMvel(COMvel);
723	>	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
724	>	}
725	>
726	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
727	>	if (states[i] == zcsMoving){
728	>	zconsMols[i]->getCOMvel(COMvel);
729	>	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
730	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
731	>	}
732	>	}
733	>
734	>	#ifndef IS_MPI
735	>	MVzOfMovingMols = MVzOfMovingMols_local;
736	>	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
737		#else
738	<	totalMZOfZCons = mzOfZCons;
738	>	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
739	>	MPI_SUM, MPI_COMM_WORLD);
740	>	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
741	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
742		#endif
743
744	<	for(int i = 0; i < unconsMols.size(); i++){
745	<	unconsMols[i]->getCOM(COM);
746	<	mzOfUncons += massOfUnconsMols[i] * COM[2];
744	>	double vzOfMovingMols;
745	>	vzOfMovingMols = MVzOfMovingMols /
746	>	(totalMassOfUncons + totalMassOfMovingZMols);
747	>
748	>	//modify the velocites of unconstrained molecules
749	>	Atom** unconsAtoms;
750	>	for (int i = 0; i < (int) (unconsMols.size()); i++){
751	>	unconsAtoms = unconsMols[i]->getMyAtoms();
752	>	for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
753	>	unconsAtoms[j]->getVel(vel);
754	>	vel[whichDirection] -= vzOfMovingMols;
755	>	unconsAtoms[j]->setVel(vel);
756	>	}
757	>	}
758	>
759	>	//modify the velocities of moving z-constrained molecuels
760	>	Atom** movingZAtoms;
761	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
762	>	if (states[i] == zcsMoving){
763	>	movingZAtoms = zconsMols[i]->getMyAtoms();
764	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
765	>	movingZAtoms[j]->getVel(vel);
766	>	vel[whichDirection] -= vzOfMovingMols;
767	>	movingZAtoms[j]->setVel(vel);
768	>	}
769	>	}
770		}
771	<
771	>
772	>
773	>	zSysCOMVel = calcSysCOMVel();
774		#ifdef IS_MPI
775	<	MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
775	>	if (worldRank == 0){
776	>	#endif
777	>	//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
778	>	#ifdef IS_MPI
779	>	}
780	>	#endif
781	>	}
782	>
783	>	/**
784	>	*
785	>	*/
786	>
787	>	template<typename T> void ZConstraint<T>::doZconstraintForce(){
788	>	Atom** zconsAtoms;
789	>	double totalFZ;
790	>	double totalFZ_local;
791	>	double COM[3];
792	>	double force[3];
793	>
794	>	//constrain the molecules which do not reach the specified positions
795	>
796	>	//Zero Out the force of z-contrained molecules
797	>	totalFZ_local = 0;
798	>
799	>	//calculate the total z-contrained force of fixed z-contrained molecules
800	>
801	>	//cout << "before zero out z-constraint force on fixed z-constraint molecuels "
802	>	//       << "total force is " << calcTotalForce() << endl;
803	>
804	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
805	>	if (states[i] == zcsFixed){
806	>	zconsMols[i]->getCOM(COM);
807	>	zconsAtoms = zconsMols[i]->getMyAtoms();
808	>
809	>	fz[i] = 0;
810	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
811	>	zconsAtoms[j]->getFrc(force);
812	>	fz[i] += force[whichDirection];
813	>	}
814	>	totalFZ_local += fz[i];
815	>
816	>	//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
817	>	//      <<"\tcurrent zpos: " << COM[whichDirection]
818	>	//      << "\tcurrent fz: " <<fz[i] << endl;
819	>	}
820	>	}
821	>
822	>	//calculate total z-constraint force
823	>	#ifdef IS_MPI
824	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
825		#else
826	<	totalMZOfUncons = mzOfUncons;
827	<	#endif
526	<
527	<	zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
826	>	totalFZ = totalFZ_local;
827	>	#endif
828
829	<	cout << "current time: " << info->getTime() <<endl;
830	<	for(int i = 0; i < zconsMols.size(); i++){
831	<
829	>
830	>	// apply negative to fixed z-constrained molecues;
831	>	force[0] = 0;
832	>	force[1] = 0;
833	>	force[2] = 0;
834	>
835	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
836	>	if (states[i] == zcsFixed){
837	>	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
838	>	zconsAtoms = zconsMols[i]->getMyAtoms();
839	>
840	>	for (int j = 0; j < nAtomOfCurZConsMol; j++){
841	>	//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
842	>	force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
843	>	zconsAtoms[j],
844	>	fz[i]);
845	>	zconsAtoms[j]->addFrc(force);
846	>	}
847	>	}
848	>	}
849	>
850	>	//cout << "after zero out z-constraint force on fixed z-constraint molecuels "
851	>	//      << "total force is " << calcTotalForce() << endl;
852	>
853	>
854	>	force[0] = 0;
855	>	force[1] = 0;
856	>	force[2] = 0;
857	>
858	>	//modify the forces of unconstrained molecules
859	>	for (int i = 0; i < (int) (unconsMols.size()); i++){
860	>	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
861	>
862	>	for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
863	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
864	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
865	>	totalFZ);
866	>	unconsAtoms[j]->addFrc(force);
867	>	}
868	>	}
869	>
870	>	//modify the forces of moving z-constrained molecules
871	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
872	>	if (states[i] == zcsMoving){
873	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
874	>
875	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
876	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
877	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
878	>	totalFZ);
879	>	movingZAtoms[j]->addFrc(force);
880	>	}
881	>	}
882	>	}
883	>	//  cout << "after substracting z-constraint force from moving molecuels "
884	>	//        << "total force is " << calcTotalForce()  << endl;
885	>	}
886	>
887	>	/**
888	>	*
889	>	*
890	>	*/
891	>
892	>	template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
893	>	double force[3];
894	>	double harmonicU;
895	>	double harmonicF;
896	>	double COM[3];
897	>	double diff;
898	>	double totalFZ_local;
899	>	double totalFZ;
900	>
901	>	force[0] = 0;
902	>	force[1] = 0;
903	>	force[2] = 0;
904	>
905	>	totalFZ_local = 0;
906	>
907	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
908	>	if (states[i] == zcsMoving){
909	>	zconsMols[i]->getCOM(COM);
910	>	//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
911	>	//     << "\tcurrent zpos: " << COM[whichDirection] << endl;
912	>
913	>	diff = COM[whichDirection] - resPos[i];
914	>
915	>	harmonicU = 0.5 * kz[i] * diff * diff;
916	>	info->lrPot += harmonicU;
917	>
918	>	harmonicF = -kz[i] * diff;
919	>	totalFZ_local += harmonicF;
920	>
921	>	//adjust force
922	>
923	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
924	>
925	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
926	>	//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
927	>	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
928	>	movingZAtoms[j],
929	>	harmonicF);
930	>	movingZAtoms[j]->addFrc(force);
931	>	}
932	>	}
933	>	}
934	>
935	>	#ifndef IS_MPI
936	>	totalFZ = totalFZ_local;
937	>	#else
938	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
939	>	#endif
940	>
941	>	//cout << "before substracting harmonic force from moving molecuels "
942	>	//      << "total force is " << calcTotalForce()  << endl;
943	>
944	>	force[0] = 0;
945	>	force[1] = 0;
946	>	force[2] = 0;
947	>
948	>	//modify the forces of unconstrained molecules
949	>	for (int i = 0; i < (int) (unconsMols.size()); i++){
950	>	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
951	>
952	>	for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
953	>	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
954	>	force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
955	>	totalFZ);
956	>	unconsAtoms[j]->addFrc(force);
957	>	}
958	>	}
959	>
960	>	//cout << "after substracting harmonic force from moving molecuels "
961	>	//      << "total force is " << calcTotalForce()  << endl;
962	>	}
963	>
964	>	/**
965	>	*
966	>	*/
967	>
968	>	template<typename T> bool ZConstraint<T>::checkZConsState(){
969	>	double COM[3];
970	>	double diff;
971	>
972	>	int changed_local;
973	>	int changed;
974	>
975	>	changed_local = 0;
976	>
977	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
978		zconsMols[i]->getCOM(COM);
979	<
980	<	cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
981	<	deltaZ = zsys + refZ[i] - COM[2];
982	<	cout << "\tdistance: " << COM[2] +deltaZ - zsys;
983	<	//update z coordinate
984	<	zconsAtoms = zconsMols[i]->getMyAtoms();
985	<	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
986	<	zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);
987	<	}
988	<
989	<	//calculate z constrain force
990	<	fz[i] = massOfZConsMols[i]* deltaZ / dt2;
991	<
992	<	cout << "\tforce: " << fz[i] << endl;
979	>	diff = fabs(COM[whichDirection] - zPos[i]);
980	>	if (diff <= zconsTol && states[i] == zcsMoving){
981	>	states[i] = zcsFixed;
982	>	changed_local = 1;
983	>
984	>	if(usingSMD)
985	>	prevCantPos = cantPos;
986	>
987	>	if (hasZConsGap)
988	>	endFixTime[i] = info->getTime() + zconsFixTime;
989	>	}
990	>	else if (diff > zconsTol && states[i] == zcsFixed){
991	>	states[i] = zcsMoving;
992	>	changed_local = 1;
993	>
994	>	if(usingSMD)
995	>	cantPos = prevCantPos;
996	>
997	>	if (hasZConsGap)
998	>	endFixTime[i] = INFINITE_TIME;
999	>	}
1000		}
1001
1002	+	#ifndef IS_MPI
1003	+	changed = changed_local;
1004	+	#else
1005	+	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1006	+	#endif
1007	+
1008	+	return (changed > 0);
1009	+	}
1010	+
1011	+	template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1012	+	int havingFixed_local;
1013	+	int havingFixed;
1014	+
1015	+	havingFixed_local = 0;
1016	+
1017	+	for (int i = 0; i < (int) (zconsMols.size()); i++)
1018	+	if (states[i] == zcsFixed){
1019	+	havingFixed_local = 1;
1020	+	break;
1021	+	}
1022	+
1023	+	#ifndef IS_MPI
1024	+	havingFixed = havingFixed_local;
1025	+	#else
1026	+	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1027	+	MPI_COMM_WORLD);
1028	+	#endif
1029	+
1030	+	return (havingFixed > 0);
1031	+	}
1032	+
1033	+
1034	+	/**
1035	+	*
1036	+	*/
1037	+	template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1038	+	int havingMoving_local;
1039	+	int havingMoving;
1040	+
1041	+	havingMoving_local = 0;
1042	+
1043	+	for (int i = 0; i < (int) (zconsMols.size()); i++)
1044	+	if (states[i] == zcsMoving){
1045	+	havingMoving_local = 1;
1046	+	break;
1047	+	}
1048	+
1049	+	#ifndef IS_MPI
1050	+	havingMoving = havingMoving_local;
1051	+	#else
1052	+	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1053	+	MPI_COMM_WORLD);
1054	+	#endif
1055	+
1056	+	return (havingMoving > 0);
1057	+	}
1058	+
1059	+	/**
1060	+	*
1061	+	*/
1062	+
1063	+	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
1064	+	double MVzOfMovingMols_local;
1065	+	double MVzOfMovingMols;
1066	+	double totalMassOfMovingZMols_local;
1067	+	double totalMassOfMovingZMols;
1068	+	double COMvel[3];
1069	+
1070	+	MVzOfMovingMols_local = 0;
1071	+	totalMassOfMovingZMols_local = 0;
1072	+
1073	+	for (int i = 0; i < unconsMols.size(); i++){
1074	+	unconsMols[i]->getCOMvel(COMvel);
1075	+	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1076	+	}
1077	+
1078	+	for (int i = 0; i < zconsMols.size(); i++){
1079	+	if (states[i] == zcsMoving){
1080	+	zconsMols[i]->getCOMvel(COMvel);
1081	+	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1082	+	totalMassOfMovingZMols_local += massOfZConsMols[i];
1083	+	}
1084	+	}
1085	+
1086	+	#ifndef IS_MPI
1087	+	MVzOfMovingMols = MVzOfMovingMols_local;
1088	+	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1089	+	#else
1090	+	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1091	+	MPI_SUM, MPI_COMM_WORLD);
1092	+	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1093	+	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1094	+	#endif
1095	+
1096	+	double vzOfMovingMols;
1097	+	vzOfMovingMols = MVzOfMovingMols /
1098	+	(totalMassOfUncons + totalMassOfMovingZMols);
1099	+
1100	+	return vzOfMovingMols;
1101	+	}
1102	+
1103	+	/**
1104	+	*
1105	+	*/
1106	+
1107	+	template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1108	+	double COMvel[3];
1109	+	double tempMVz_local;
1110	+	double tempMVz;
1111	+	double massOfZCons_local;
1112	+	double massOfZCons;
1113	+
1114	+
1115	+	tempMVz_local = 0;
1116	+
1117	+	for (int i = 0 ; i < nMols; i++){
1118	+	molecules[i].getCOMvel(COMvel);
1119	+	tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
1120	+	}
1121	+
1122	+	massOfZCons_local = 0;
1123	+
1124	+	for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
1125	+	massOfZCons_local += massOfZConsMols[i];
1126	+	}
1127	+	#ifndef IS_MPI
1128	+	massOfZCons = massOfZCons_local;
1129	+	tempMVz = tempMVz_local;
1130	+	#else
1131	+	MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
1132	+	MPI_COMM_WORLD);
1133	+	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1134	+	#endif
1135	+
1136	+	return tempMVz / (totalMassOfUncons + massOfZCons);
1137	+	}
1138	+
1139	+	/**
1140	+	*
1141	+	*/
1142	+
1143	+	template<typename T> double ZConstraint<T>::calcTotalForce(){
1144	+	double force[3];
1145	+	double totalForce_local;
1146	+	double totalForce;
1147	+
1148	+	totalForce_local = 0;
1149	+
1150	+	for (int i = 0; i < nAtoms; i++){
1151	+	atoms[i]->getFrc(force);
1152	+	totalForce_local += force[whichDirection];
1153	+	}
1154	+
1155	+	#ifndef IS_MPI
1156	+	totalForce = totalForce_local;
1157	+	#else
1158	+	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1159	+	MPI_COMM_WORLD);
1160	+	#endif
1161	+
1162	+	return totalForce;
1163	+	}
1164	+
1165	+	/**
1166	+	*
1167	+	*/
1168	+
1169	+	template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1170	+	//calculate the number of atoms of moving z-constrained molecules
1171	+	int nMovingZAtoms_local;
1172	+	int nMovingZAtoms;
1173	+
1174	+	nMovingZAtoms_local = 0;
1175	+	for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1176	+	if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1177	+	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1178	+	}
1179	+
1180	+	#ifdef IS_MPI
1181	+	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
1182	+	MPI_COMM_WORLD);
1183	+	#else
1184	+	nMovingZAtoms = nMovingZAtoms_local;
1185	+	#endif
1186	+	totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1187	+	}
1188	+
1189	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
1190	+	Atom* atom,
1191	+	double totalForce){
1192	+	return totalForce / mol->getNAtoms();
1193	+	}
1194	+
1195	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
1196	+	double totalForce){
1197	+	return totalForce / totNumOfMovingAtoms;
1198	+	}
1199	+
1200	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
1201	+	Atom* atom,
1202	+	double totalForce){
1203	+	return totalForce / mol->getNAtoms();
1204	+	}
1205	+
1206	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1207	+	double totalForce){
1208	+	return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1209	+	}
1210	+
1211	+	/**
1212	+	*
1213	+	*/
1214	+
1215	+	template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1216	+	//calculate the number of atoms of moving z-constrained molecules
1217	+	double massOfMovingZAtoms_local;
1218	+	double massOfMovingZAtoms;
1219	+
1220	+	massOfMovingZAtoms_local = 0;
1221	+	for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1222	+	if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1223	+	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1224	+	}
1225	+
1226	+	#ifdef IS_MPI
1227	+	MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
1228	+	MPI_SUM, MPI_COMM_WORLD);
1229	+	#else
1230	+	massOfMovingZAtoms = massOfMovingZAtoms_local;
1231	+	#endif
1232	+	totMassOfMovingAtoms = massOfMovingZAtoms +
1233	+	zconsIntegrator->totalMassOfUncons;
1234	+	}
1235	+
1236	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
1237	+	Atom* atom,
1238	+	double totalForce){
1239	+	return totalForce * atom->getMass() / mol->getTotalMass();
1240	+	}
1241	+
1242	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
1243	+	double totalForce){
1244	+	return totalForce * atom->getMass() / totMassOfMovingAtoms;
1245	+	}
1246	+
1247	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
1248	+	Atom* atom,
1249	+	double totalForce){
1250	+	return totalForce * atom->getMass() / mol->getTotalMass();
1251	+	}
1252	+
1253	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
1254	+	double totalForce){
1255	+	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1256	+	}
1257	+
1258	+	template<typename T> void ZConstraint<T>::updateZPos(){
1259	+	double curTime;
1260	+	double COM[3];
1261	+
1262	+	curTime = info->getTime();
1263	+
1264	+	for (size_t i = 0; i < zconsMols.size(); i++){
1265	+
1266	+	if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1267	+	zPos[i] += zconsGap;
1268	+
1269	+	if (usingSMD){
1270	+	zconsMols[i]->getCOM(COM);
1271	+	cantPos[i] = COM[whichDirection];
1272	+	}
1273
1274	+	}
1275	+
1276	+	}
1277	+
1278		}
1279	+
1280	+	template<typename T> void ZConstraint<T>::updateCantPos(){
1281	+	double curTime;
1282	+	double dt;
1283	+
1284	+	curTime = info->getTime();
1285	+	dt = info->dt;
1286	+
1287	+	for (size_t i = 0; i < zconsMols.size(); i++){
1288	+	if (states[i] == zcsMoving){
1289	+	cantPos[i] += cantVel[i] * dt;
1290	+	}
1291	+	}
1292	+
1293	+	}

Diff Legend

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