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 829 by gezelter, Tue Oct 28 16:03:37 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	–	template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5	–	: T(theInfo, the_ff), indexOfZConsMols(NULL), fz(NULL), curZPos(NULL),
6	–	fzOut(NULL), curZconsTime(0), forcePolicy(NULL)
7	–	{
4
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		ZConsParaData* zConsParaData;
16		DoubleData* sampleTime;
17		DoubleData* tolerance;
18	+	DoubleData* gap;
19	+	DoubleData* fixtime;
20		StringData* policy;
21		StringData* filename;
22	+	IntData* smdFlag;
23		double COM[3];
24
25		//by default, the direction of constraint is z
#	Line 23 | Line 30 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
30
31		//estimate the force constant of harmonical potential
32		double Kb = 1.986E-3 ; //in kcal/K
26	–
27	–	double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
28	–	zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
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");
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);
47	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
48		}
49		else{
50		policy = dynamic_cast<StringData*>(data);
51	<
52	<	if(!policy){
53	<	sprintf( painCave.errMsg,
54	<	"ZConstraint Error: Convertion from GenericData to StringData failure, "
55	<	"PolicyByMass is used\n");
51	>
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);
59	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
60		}
61		else{
62	<	if(policy->getData() == "BYNUMBER")
63	<	forcePolicy = (ForceSubtractionPolicy*) new PolicyByNumber(this);
64	<	else if(policy->getData() == "BYMASS")
65	<	forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this);
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");
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	<	}
72	>	forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
73	>	}
74		}
75		}
76	<
77	<
76	>
77	>
78		//retrieve sample time of z-contraint
79		data = info->getProperty(ZCONSTIME_ID);
80	<
81	<	if(!data) {
82	<
83	<	sprintf( painCave.errMsg,
84	<	"ZConstraint error: If you use an ZConstraint\n"
77	<	" , you must set sample time.\n");
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{
82	–
89		sampleTime = dynamic_cast<DoubleData*>(data);
84	–
85	–	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();
91	<
95	>	simError();
96		}
97		else{
98		this->zconsTime = sampleTime->getData();
99		}
96	–
100		}
101	<
101	>
102		//retrieve output filename of z force
103		data = info->getProperty(ZCONSFILENAME_ID);
104	<	if(!data) {
105	<
106	<
107	<	sprintf( painCave.errMsg,
105	<	"ZConstraint error: If you use an ZConstraint\n"
106	<	" , you must set output filename of z-force.\n");
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();
109	<
109	>	simError();
110		}
111		else{
112	–
112		filename = dynamic_cast<StringData*>(data);
114	–
115	–	if(!filename){
113
114	<	sprintf( painCave.errMsg,
115	<	"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();
121	<
118	>	simError();
119		}
120		else{
121		this->zconsOutput = filename->getData();
122		}
126	–
123		}
124
125		//retrieve tolerance for z-constraint molecuels
126		data = info->getProperty(ZCONSTOL_ID);
127	<
128	<	if(!data) {
129	<
134	<	sprintf( painCave.errMsg,
135	<	"ZConstraint error: can not get tolerance \n");
127	>
128	>	if (!data){
129	>	sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n");
130		painCave.isFatal = 1;
131	<	simError();
131	>	simError();
132		}
133		else{
140	–
134		tolerance = dynamic_cast<DoubleData*>(data);
142	–
143	–	if(!tolerance){
135
136	<	sprintf( painCave.errMsg,
137	<	"ZConstraint error: Can not get property from SimInfo\n");
136	>	if (!tolerance){
137	>	sprintf(painCave.errMsg,
138	>	"ZConstraint error: Can not get property from SimInfo\n");
139		painCave.isFatal = 1;
140	<	simError();
149	<
140	>	simError();
141		}
142		else{
143		this->zconsTol = tolerance->getData();
144		}
145	+	}
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		}
156	–
157	–	//retrieve index of z-constraint molecules
158	–	data = info->getProperty(ZCONSPARADATA_ID);
159	–	if(!data) {
164
165	<	sprintf( painCave.errMsg,
166	<	"ZConstraint error: If you use an ZConstraint\n"
167	<	" , you must set index of z-constraint molecules.\n");
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	>	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	>
188	>
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	>
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();
216	>	simError();
217		}
218		else{
168	–
219		zConsParaData = dynamic_cast<ZConsParaData*>(data);
170	–
171	–	if(!zConsParaData){
220
221	<	sprintf( painCave.errMsg,
222	<	"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
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();
177	<
225	>	simError();
226		}
227		else{
180	–
228		parameters = zConsParaData->getData();
229
230		//check the range of zconsIndex
#	Line 189 | Line 236 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
236		int totalNumMol;
237
238		minIndex = (*parameters)[0].zconsIndex;
239	<	if(minIndex < 0){
240	<	sprintf( painCave.errMsg,
194	<	"ZConstraint error: index is out of range\n");
239	>	if (minIndex < 0){
240	>	sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
241		painCave.isFatal = 1;
242	<	simError();
243	<	}
242	>	simError();
243	>	}
244
245		maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
246
247		#ifndef IS_MPI
248		totalNumMol = nMols;
249		#else
250	<	totalNumMol = mpiSim->getTotNmol();
250	>	totalNumMol = mpiSim->getNMolGlobal();
251		#endif
252	<
253	<	if(maxIndex > totalNumMol - 1){
254	<	sprintf( painCave.errMsg,
209	<	"ZConstraint error: index is out of range\n");
252	>
253	>	if (maxIndex > totalNumMol - 1){
254	>	sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n");
255		painCave.isFatal = 1;
256	<	simError();
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	<
218	<	if(!(*parameters)[i].havingZPos){
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++){
264	>	for (int j = 0; j < nMols; j++){
265		if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
266		molecules[j].getCOM(COM);
267	<	break;
267	>	break;
268		}
269		}
270		#else
271	<	//query which processor current zconstraint molecule belongs to
272	<	int *MolToProcMap;
271	>	//query which processor current zconstraint molecule belongs to
272	>	int* MolToProcMap;
273		int whichNode;
274	<
274	>
275		MolToProcMap = mpiSim->getMolToProcMap();
276		whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
277	<
277	>
278		//broadcast the zpos of current z-contraint molecule
279		//the node which contain this
280	<
281	<	if (worldRank == whichNode ){
282	<
239	<	for(int j = 0; j < nMols; j++)
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;
285	>	break;
286		}
244	–
287		}
288
289	<	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);
289	>	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode,
290	>	MPI_COMM_WORLD);
291		#endif
292	<
292	>
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	<
258	<	}
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		}
260	–
302		}//end if (!zConsParaData)
303	+
304		}//end  if (!data)
305	<
306	<	//
305	>
306	>	//
307		#ifdef IS_MPI
308		update();
309		#else
310		int searchResult;
311	<
312	<	for(int i = 0; i < nMols; i++){
271	<
311	>
312	>	for (int i = 0; i < nMols; i++){
313		searchResult = isZConstraintMol(&molecules[i]);
314	<
315	<	if(searchResult > -1){
275	<
314	>
315	>	if (searchResult > -1){
316		zconsMols.push_back(&molecules[i]);
317		massOfZConsMols.push_back(molecules[i].getTotalMass());
318
319		zPos.push_back((*parameters)[searchResult].zPos);
280	–	//       cout << "index: "<< (*parameters)[searchResult].zconsIndex
281	–	//              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
282	–
320		kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
321	<	molecules[i].getCOM(COM);
321	>
322	>	if(usingSMD)
323	>	cantVel.push_back((*parameters)[searchResult].cantVel);
324	>
325		}
326	<	else
287	<	{
288	<
326	>	else{
327		unconsMols.push_back(&molecules[i]);
328		massOfUnconsMols.push_back(molecules[i].getTotalMass());
291	–
329		}
330		}
331
332	<	fz = new double[zconsMols.size()];
333	<	curZPos = new double[zconsMols.size()];
334	<	indexOfZConsMols = new int [zconsMols.size()];
332	>	fz.resize(zconsMols.size());
333	>	curZPos.resize(zconsMols.size());
334	>	indexOfZConsMols.resize(zconsMols.size());
335
299	–	if(!fz || !curZPos || !indexOfZConsMols){
300	–	sprintf( painCave.errMsg,
301	–	"Memory allocation failure in class Zconstraint\n");
302	–	painCave.isFatal = 1;
303	–	simError();
304	–	}
305	–
336		//determine the states of z-constraint molecules
337	<	for(int i = 0; i < (int)(zconsMols.size()); i++){
337	>	for (size_t i = 0; i < zconsMols.size(); i++){
338		indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
339
340		zconsMols[i]->getCOM(COM);
341	<	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
341	>
342	>	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
343		states.push_back(zcsFixed);
344	<	else
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	<
358	>
359	>	if(usingSMD)
360	>	prevCantPos = cantPos;
361		#endif
362
363	+
364		//get total masss of unconstraint molecules
365		double totalMassOfUncons_local;
366		totalMassOfUncons_local = 0;
367	<
368	<	for(int i = 0; i < (int)(unconsMols.size()); i++)
367	>
368	>	for (size_t i = 0; i < unconsMols.size(); i++)
369		totalMassOfUncons_local += unconsMols[i]->getTotalMass();
370	<
370	>
371		#ifndef IS_MPI
372		totalMassOfUncons = totalMassOfUncons_local;
373		#else
374	<	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
375	<	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
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++)
381	>	for (int i = 0; i < (int) (unconsMols.size()); i++)
382		nUnconsAtoms_local += unconsMols[i]->getNAtoms();
383	<
383	>
384		#ifndef IS_MPI
385		totNumOfUnconsAtoms = nUnconsAtoms_local;
386		#else
387	<	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1,
388	<	MPI_INT,MPI_SUM, MPI_COMM_WORLD);
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()
350	<	{
351	<	if(fz)
352	<	delete[] fz;
394	>	template<typename T> ZConstraint<T>::~ZConstraint(){
395
396	<	if(curZPos)
355	<	delete[] curZPos;
356	<
357	<	if(indexOfZConsMols)
358	<	delete[] indexOfZConsMols;
359	<
360	<	if(fzOut)
396	>	if (fzOut){
397		delete fzOut;
398	<
399	<	if(forcePolicy)
398	>	}
399	>
400	>	if (forcePolicy){
401		delete forcePolicy;
402	+	}
403		}
404
405
#	Line 370 | Line 408 | template<typename T> void ZConstraint<T>::update()
408		*/
409
410		#ifdef IS_MPI
411	<	template<typename T> void ZConstraint<T>::update()
374	<	{
411	>	template<typename T> void ZConstraint<T>::update(){
412		double COM[3];
413		int index;
414	<
414	>
415		zconsMols.clear();
416		massOfZConsMols.clear();
417		zPos.clear();
418		kz.clear();
419	<
419	>	cantPos.clear();
420	>	cantVel.clear();
421	>
422		unconsMols.clear();
423		massOfUnconsMols.clear();
385	–
424
425	+
426		//creat zconsMol and unconsMol lists
427	<	for(int i = 0; i < nMols; i++){
389	<
427	>	for (int i = 0; i < nMols; i++){
428		index = isZConstraintMol(&molecules[i]);
429	<
430	<	if(index > -1){
393	<
429	>
430	>	if (index > -1){
431		zconsMols.push_back(&molecules[i]);
432		zPos.push_back((*parameters)[index].zPos);
433		kz.push_back((*parameters)[index].kRatio * zForceConst);
434	<	massOfZConsMols.push_back(molecules[i].getTotalMass());
435	<
436	<	molecules[i].getCOM(COM);
434	>	massOfZConsMols.push_back(molecules[i].getTotalMass());
435	>
436	>	if(usingSMD)
437	>	cantVel.push_back((*parameters)[index].cantVel);
438	>
439		}
440	<	else
402	<	{
403	<
440	>	else{
441		unconsMols.push_back(&molecules[i]);
442		massOfUnconsMols.push_back(molecules[i].getTotalMass());
406	–
443		}
444		}
445
446	<	//determine the states of z-constraint molecules
447	<	for(int i = 0; i < (int)(zconsMols.size()); i++){
448	<	zconsMols[i]->getCOM(COM);
449	<	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
450	<	states.push_back(zcsFixed);
451	<	else
416	<	states.push_back(zcsMoving);
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	+	//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	<	//The reason to declare fz and indexOfZconsMols as pointer to array is
460	<	// that we want to make the MPI communication simple
461	<	if(fz)
462	<	delete[] fz;
463	<
425	<	if(curZPos)
426	<	delete[] curZPos;
427	<
428	<	if(indexOfZConsMols)
429	<	delete[] indexOfZConsMols;
430	<
431	<	if (zconsMols.size() > 0){
432	<	fz = new double[zconsMols.size()];
433	<	curZPos = new double[zconsMols.size()];
434	<	indexOfZConsMols =  new int[zconsMols.size()];
435	<
436	<	if(!fz || !curZPos || !indexOfZConsMols){
437	<	sprintf( painCave.errMsg,
438	<	"Memory allocation failure in class Zconstraint\n");
439	<	painCave.isFatal = 1;
440	<	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 < (int)(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	<	curZPos = NULL;
476	<	indexOfZConsMols = NULL;
477	<	}
478	<
472	>	if(usingSMD)
473	>	cantPos.push_back(COM[whichDirection]);
474	>	}
475	>
476	>	if(usingSMD)
477	>	prevCantPos = cantPos;
478	>
479		//
480		forcePolicy->update();
456	–
481		}
482
483		#endif
#	Line 467 | Line 491 | template<typename T> int ZConstraint<T>::isZConstraint
491		*    other non-negative values, its index in indexOfAllZConsMols vector
492		*/
493
494	<	template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
471	<	{
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 = parameters->size() - 1;
504	<
504	>
505		//Binary Search (we have sorted the array)
506	<	while(low <= high){
507	<	mid = (low + high) /2;
506	>	while (low <= high){
507	>	mid = (low + high) / 2;
508		if ((*parameters)[mid].zconsIndex == index)
509		return mid;
510	<	else if ((*parameters)[mid].zconsIndex > 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		template<typename T> void ZConstraint<T>::integrate(){
497	–
520		// creat zconsWriter
521		fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
522	<
523	<	if(!fzOut){
524	<	sprintf( painCave.errMsg,
503	<	"Memory allocation failure in class Zconstraint\n");
522	>
523	>	if (!fzOut){
524	>	sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
525		painCave.isFatal = 1;
526		simError();
527		}
528	<
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();
513	–
514	–	T::integrate();
534
535	+	T::integrate();
536		}
517	–
537
538	+
539		/**
540		*
541		*
#	Line 530 | Line 550 | template<typename T> void ZConstraint<T>::calcForce(in
550
551		T::calcForce(calcPot, calcStress);
552
553	<	if (checkZConsState()){
553	>
554	>	if (hasZConsGap){
555	>	updateZPos();
556	>	}
557	>
558	>	if (checkZConsState()){
559		zeroOutVel();
560		forcePolicy->update();
561		}
562	<
562	>
563		zsys = calcZSys();
564		zSysCOMVel = calcSysCOMVel();
565		#ifdef IS_MPI
566	<	if(worldRank == 0){
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;
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())
578	>	if (haveFixedZMols()){
579		this->doZconstraintForce();
580	+	}
581
582	<	//use harmonical poteintial to move the molecules to the specified positions
583	<	if (haveMovingZMols())
584	<	this->doHarmonic();
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++){
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];
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];
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	<	}
608	<	fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
609	<	curZconsTime += zconsTime;
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){
614	>	if (worldRank == 0){
615		#endif
616	<	//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
616	>	//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
617		#ifdef IS_MPI
618		}
619		#endif
#	Line 592 | Line 623 | template<typename T> void ZConstraint<T>::calcForce(in
623		/**
624		*
625		*/
626	<
627	<	template<typename T> double ZConstraint<T>::calcZSys()
597	<	{
626	>
627	>	template<typename T> double ZConstraint<T>::calcZSys(){
628		//calculate reference z coordinate for z-constraint molecules
629		double totalMass_local;
630		double totalMass;
#	Line 602 | Line 632 | template<typename T> double ZConstraint<T>::calcZSys()
632		double totalMZ;
633		double massOfCurMol;
634		double COM[3];
635	<
635	>
636		totalMass_local = 0;
637		totalMZ_local = 0;
638	<
639	<	for(int i = 0; i < nMols; i++){
638	>
639	>	for (int i = 0; i < nMols; i++){
640		massOfCurMol = molecules[i].getTotalMass();
641		molecules[i].getCOM(COM);
642	<
642	>
643		totalMass_local += massOfCurMol;
644		totalMZ_local += massOfCurMol * COM[whichDirection];
615	–
645		}
646	<
647	<
646	>
647	>
648		#ifdef IS_MPI
649	<	MPI_Allreduce(&totalMass_local, &totalMass, 1,
650	<	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
651	<	MPI_Allreduce(&totalMZ_local, &totalMZ, 1,
652	<	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
624	<	#else
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
655	>	#endif
656
657		double zsys;
658		zsys = totalMZ / totalMass;
#	Line 635 | Line 663 | template<typename T> void ZConstraint<T>::thermalize(
663		/**
664		*
665		*/
666	<	template<typename T> void ZConstraint<T>::thermalize( void ){
639	<
666	>	template<typename T> void ZConstraint<T>::thermalize(void){
667		T::thermalize();
668		zeroOutVel();
669		}
#	Line 646 | Line 673 | template<typename T> void ZConstraint<T>::zeroOutVel()
673		*/
674
675		template<typename T> void ZConstraint<T>::zeroOutVel(){
649	–
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
656	–
657	–	for(int i = 0; i < (int)(zconsMols.size()); i++){
682
683	<	if (states[i] == zcsFixed){
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
661	–	zconsMols[i]->getCOMvel(COMvel);
662	–	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
663	–
688		fixedZAtoms = zconsMols[i]->getMyAtoms();
689	<
690	<	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
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);
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;
696	>	zconsMols[i]->getCOMvel(COMvel);
697	>	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
698		}
675	–
699		}
700
701	<	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
701	>	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
702
703		zSysCOMVel = calcSysCOMVel();
704		#ifdef IS_MPI
705	<	if(worldRank == 0){
705	>	if (worldRank == 0){
706		#endif
707	<	//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
707	>	//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
708		#ifdef IS_MPI
709		}
710		#endif
711	<
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	<
717	>
718		MVzOfMovingMols_local = 0;
719		totalMassOfMovingZMols_local = 0;
720
721	<	for(int i =0; i < (int)(unconsMols.size()); i++){
721	>	for (int i = 0; i < (int) (unconsMols.size()); i++){
722		unconsMols[i]->getCOMvel(COMvel);
723	<	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
723	>	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
724		}
725
726	<	for(int i = 0; i < (int)(zconsMols.size()); i++){
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];
730	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
731		}
709	–
732		}
733
734		#ifndef IS_MPI
735		MVzOfMovingMols = MVzOfMovingMols_local;
736		totalMassOfMovingZMols = totalMassOfMovingZMols_local;
737		#else
738	<	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
739	<	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
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		double vzOfMovingMols;
745	<	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
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++){
726	<
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++){
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		}
733	–
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	<
740	<	if (states[i] ==zcsMoving){
741	<
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++){
764	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
765		movingZAtoms[j]->getVel(vel);
766		vel[whichDirection] -= vzOfMovingMols;
767	<	movingZAtoms[j]->setVel(vel);
767	>	movingZAtoms[j]->setVel(vel);
768	>	}
769		}
770	<
749	<	}
770	>	}
771
751	–	}
772
753	–
773		zSysCOMVel = calcSysCOMVel();
774		#ifdef IS_MPI
775	<	if(worldRank == 0){
775	>	if (worldRank == 0){
776		#endif
777	<	//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
777	>	//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
778		#ifdef IS_MPI
779		}
780		#endif
762	–
781		}
782
783		/**
#	Line 767 | Line 785 | template<typename T> void ZConstraint<T>::doZconstrain
785		*/
786
787		template<typename T> void ZConstraint<T>::doZconstraintForce(){
770	–
788		Atom** zconsAtoms;
789		double totalFZ;
790		double totalFZ_local;
#	Line 775 | Line 792 | template<typename T> void ZConstraint<T>::doZconstrain
792		double force[3];
793
794		//constrain the molecules which do not reach the specified positions
795	<
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	<
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++){
788	<
804	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
805		if (states[i] == zcsFixed){
790	–
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++) {
810	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
811		zconsAtoms[j]->getFrc(force);
812	<	fz[i] += force[whichDirection];
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;
804	–
805	–
819		}
807	–
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);
824	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
825		#else
826		totalFZ = totalFZ_local;
827		#endif
828
829	<
829	>
830		// apply negative to fixed z-constrained molecues;
831	<	force[0]= 0;
832	<	force[1]= 0;
833	<	force[2]= 0;
831	>	force[0] = 0;
832	>	force[1] = 0;
833	>	force[2] = 0;
834
835	<	for(int i = 0; i < (int)(zconsMols.size()); i++){
836	<
825	<	if (states[i] == zcsFixed){
826	<
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++) {
839	>
840	>	for (int j = 0; j < nAtomOfCurZConsMol; j++){
841		//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
842	<	force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
842	>	force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
843	>	zconsAtoms[j],
844	>	fz[i]);
845		zconsAtoms[j]->addFrc(force);
846		}
835	–
847		}
837	–
848		}
849
850		//cout << "after zero out z-constraint force on fixed z-constraint molecuels "
851		//      << "total force is " << calcTotalForce() << endl;
842	–
852
844	–	force[0]= 0;
845	–	force[1]= 0;
846	–	force[2]= 0;
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	<
861	<	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
862	<
863	<	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
864	<	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
865	<	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
866	<	unconsAtoms[j]->addFrc(force);
867	<	}
858	<
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++) {
870	>	//modify the forces of moving z-constrained molecules
871	>	for (int i = 0; i < (int) (zconsMols.size()); i++){
872		if (states[i] == zcsMoving){
864	–
873		Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
874
875	<	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
875	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
876		//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
877	<	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
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;
876	<
883	>	//  cout << "after substracting z-constraint force from moving molecuels "
884	>	//        << "total force is " << calcTotalForce()  << endl;
885		}
886
887		/**
#	Line 881 | Line 889 | template<typename T> void ZConstraint<T>::doHarmonic()
889		*
890		*/
891
892	<	template<typename T> void ZConstraint<T>::doHarmonic(){
892	>	template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
893		double force[3];
894		double harmonicU;
895		double harmonicF;
#	Line 889 | Line 897 | template<typename T> void ZConstraint<T>::doHarmonic()
897		double diff;
898		double totalFZ_local;
899		double totalFZ;
900	<
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++) {
900	<
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;
910	>	//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
911	>	//     << "\tcurrent zpos: " << COM[whichDirection] << endl;
912
913	<	diff = COM[whichDirection] -zPos[i];
914	<
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;
918	>	harmonicF = -kz[i] * diff;
919		totalFZ_local += harmonicF;
920
921	<	//adjust force
922	<
921	>	//adjust force
922	>
923		Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
924
925	<	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
925	>	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
926		//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
927	<	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
928	<	movingZAtoms[j]->addFrc(force);
929	<	}
927	>	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
928	>	movingZAtoms[j],
929	>	harmonicF);
930	>	movingZAtoms[j]->addFrc(force);
931	>	}
932		}
924	–
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);
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;
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	<
951	<	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
952	<
953	<	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
954	<	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
955	<	force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
956	<	unconsAtoms[j]->addFrc(force);
957	<	}
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;
954	–
962		}
963
964		/**
#	Line 961 | Line 968 | template<typename T> bool ZConstraint<T>::checkZConsSt
968		template<typename T> bool ZConstraint<T>::checkZConsState(){
969		double COM[3];
970		double diff;
971	<
971	>
972		int changed_local;
973		int changed;
974	<
974	>
975		changed_local = 0;
969	–
970	–	for(int i =0; i < (int)(zconsMols.size()); i++){
976
977	+	for (int i = 0; i < (int) (zconsMols.size()); i++){
978		zconsMols[i]->getCOM(COM);
979		diff = fabs(COM[whichDirection] - zPos[i]);
980	<	if (  diff <= zconsTol && states[i] == zcsMoving){
980	>	if (diff <= zconsTol && states[i] == zcsMoving){
981		states[i] = zcsFixed;
982	<	changed_local = 1;
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){
990	>	else if (diff > zconsTol && states[i] == zcsFixed){
991		states[i] = zcsMoving;
992	<	changed_local = 1;
992	>	changed_local = 1;
993	>
994	>	if(usingSMD)
995	>	cantPos = prevCantPos;
996	>
997	>	if (hasZConsGap)
998	>	endFixTime[i] = INFINITE_TIME;
999		}
982	–
1000		}
1001
1002		#ifndef IS_MPI
1003	<	changed =changed_local;
1003	>	changed = changed_local;
1004		#else
1005	<	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1005	>	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1006		#endif
1007
1008		return (changed > 0);
992	–
1009		}
1010
1011		template<typename T> bool ZConstraint<T>::haveFixedZMols(){
996	–
1012		int havingFixed_local;
1013		int havingFixed;
1014
1015		havingFixed_local = 0;
1016
1017	<	for(int i = 0; i < (int)(zconsMols.size()); i++)
1017	>	for (int i = 0; i < (int) (zconsMols.size()); i++)
1018		if (states[i] == zcsFixed){
1019		havingFixed_local = 1;
1020	<	break;
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, MPI_COMM_WORLD);
1026	>	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1027	>	MPI_COMM_WORLD);
1028		#endif
1029
1030		return (havingFixed > 0);
#	Line 1019 | Line 1035 | template<typename T> bool ZConstraint<T>::haveMovingZM
1035		*
1036		*/
1037		template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1022	–
1038		int havingMoving_local;
1039		int havingMoving;
1040
1041		havingMoving_local = 0;
1042
1043	<	for(int i = 0; i < (int)(zconsMols.size()); i++)
1043	>	for (int i = 0; i < (int) (zconsMols.size()); i++)
1044		if (states[i] == zcsMoving){
1045		havingMoving_local = 1;
1046	<	break;
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, MPI_COMM_WORLD);
1052	>	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1053	>	MPI_COMM_WORLD);
1054		#endif
1055
1056		return (havingMoving > 0);
1041	–
1057		}
1058
1059		/**
1060		*
1061		*/
1062
1063	<	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1049	<	{
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	<
1069	>
1070		MVzOfMovingMols_local = 0;
1071		totalMassOfMovingZMols_local = 0;
1072
1073	<	for(int i =0; i < unconsMols.size(); i++){
1073	>	for (int i = 0; i < unconsMols.size(); i++){
1074		unconsMols[i]->getCOMvel(COMvel);
1075	<	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1075	>	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1076		}
1077
1078	<	for(int i = 0; i < zconsMols.size(); i++){
1065	<
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];
1082	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
1083		}
1071	–
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,MPI_SUM, MPI_COMM_WORLD);
1091	<	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
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 / (totalMassOfUncons + totalMassOfMovingZMols);
1097	>	vzOfMovingMols = MVzOfMovingMols /
1098	>	(totalMassOfUncons + totalMassOfMovingZMols);
1099
1100		return vzOfMovingMols;
1101		}
#	Line 1089 | Line 1104 | template<typename T> double ZConstraint<T>::calcSysCOM
1104		*
1105		*/
1106
1107	<	template<typename T> double ZConstraint<T>::calcSysCOMVel()
1093	<	{
1107	>	template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1108		double COMvel[3];
1109		double tempMVz_local;
1110		double tempMVz;
#	Line 1098 | Line 1112 | template<typename T> double ZConstraint<T>::calcSysCOM
1112		double massOfZCons;
1113
1114
1115	<	tempMVz_local = 0;
1115	>	tempMVz_local = 0;
1116
1117	<	for(int i =0 ; i < nMols; i++){
1117	>	for (int i = 0 ; i < nMols; i++){
1118		molecules[i].getCOMvel(COMvel);
1119	<	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
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++){
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, MPI_COMM_WORLD);
1132	<	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
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);
1136	>	return tempMVz / (totalMassOfUncons + massOfZCons);
1137		}
1138
1139		/**
#	Line 1126 | Line 1141 | template<typename T> double ZConstraint<T>::calcTotalF
1141		*/
1142
1143		template<typename T> double ZConstraint<T>::calcTotalForce(){
1129	–
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++){
1150	>	for (int i = 0; i < nAtoms; i++){
1151		atoms[i]->getFrc(force);
1152		totalForce_local += force[whichDirection];
1153		}
#	Line 1141 | Line 1155 | template<typename T> double ZConstraint<T>::calcTotalF
1155		#ifndef IS_MPI
1156		totalForce = totalForce_local;
1157		#else
1158	<	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1158	>	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1159	>	MPI_COMM_WORLD);
1160		#endif
1161
1162		return totalForce;
1148	–
1163		}
1164
1165		/**
#	Line 1156 | Line 1170 | template<typename T> void ZConstraint<T>::PolicyByNumb
1170		//calculate the number of atoms of moving z-constrained molecules
1171		int nMovingZAtoms_local;
1172		int nMovingZAtoms;
1173	<
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	<
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, MPI_COMM_WORLD);
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, Atom* atom, double totalForce){
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, double totalForce){
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, Atom* atom, double totalForce){
1201	<	return totalForce / mol->getNAtoms();
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, double totalForce){
1206	>	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1207	>	double totalForce){
1208		return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1209		}
1210
#	Line 1194 | Line 1216 | template<typename T> void ZConstraint<T>::PolicyByMass
1216		//calculate the number of atoms of moving z-constrained molecules
1217		double massOfMovingZAtoms_local;
1218		double massOfMovingZAtoms;
1219	<
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	<
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,MPI_SUM, MPI_COMM_WORLD);
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 + zconsIntegrator->totalMassOfUncons;
1232	>	totMassOfMovingAtoms = massOfMovingZAtoms +
1233	>	zconsIntegrator->totalMassOfUncons;
1234		}
1235
1236	<	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
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, double totalForce){
1243	<	return totalForce * atom->getMass() / totMassOfMovingAtoms;
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, Atom* atom, double totalForce){
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, double totalForce){
1254	<	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
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