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 1093 by tim, Wed Mar 17 14:22:59 2004 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	<	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
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	<	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");
56	<	painCave.isFatal = 0;
57	<	simError();
58	<
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);
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	<	}
74	<	}
75	<	}
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();
87	<	}
88	<	else{
89	<	sampleTime = dynamic_cast<DoubleData*>(data);
90	<
91	<	if (!sampleTime){
92	<	sprintf(painCave.errMsg,
93	<	"ZConstraint error: Can not get property from SimInfo\n");
94	<	painCave.isFatal = 1;
95	<	simError();
96	<	}
97	<	else{
98	<	this->zconsTime = sampleTime->getData();
99	<	}
100	<	}
101	<
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();
110	<	}
111	<	else{
112	<	filename = dynamic_cast<StringData*>(data);
113	<
114	<	if (!filename){
115	<	sprintf(painCave.errMsg,
116	<	"ZConstraint error: Can not get property from SimInfo\n");
117	<	painCave.isFatal = 1;
118	<	simError();
119	<	}
120	<	else{
121	<	this->zconsOutput = filename->getData();
122	<	}
123	<	}
124	<
125	<	//retrieve tolerance for z-constraint molecuels
126	<	data = info->getProperty(ZCONSTOL_ID);
127	<
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	<	if (!tolerance){
137	<	sprintf(painCave.errMsg,
138	<	"ZConstraint error: Can not get property from SimInfo\n");
139	<	painCave.isFatal = 1;
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	<	}
164	<
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();
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->getTotNmol();
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	<	(*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;
311	<
312	<	for (int i = 0; i < nMols; i++){
313	<	searchResult = isZConstraintMol(&molecules[i]);
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);
320	<	kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
321	<
322	<	if(usingSMD)
323	<	cantVel.push_back((*parameters)[searchResult].cantVel);
324	<
325	<	}
326	<	else{
327	<	unconsMols.push_back(&molecules[i]);
328	<	massOfUnconsMols.push_back(molecules[i].getTotalMass());
329	<	}
330	<	}
331	<
332	<	fz.resize(zconsMols.size());
333	<	curZPos.resize(zconsMols.size());
334	<	indexOfZConsMols.resize(zconsMols.size());
335	<
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	<	#endif
360	<
361	<
362	<	//get total masss of unconstraint molecules
363	<	double totalMassOfUncons_local;
364	<	totalMassOfUncons_local = 0;
365	<
366	<	for (size_t i = 0; i < unconsMols.size(); i++)
367	<	totalMassOfUncons_local += unconsMols[i]->getTotalMass();
368	<
369	<	#ifndef IS_MPI
370	<	totalMassOfUncons = totalMassOfUncons_local;
371	<	#else
372	<	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,
373	<	MPI_SUM, MPI_COMM_WORLD);
374	<	#endif
375	<
376	<	//get total number of unconstrained atoms
377	<	int nUnconsAtoms_local;
378	<	nUnconsAtoms_local = 0;
379	<	for (int i = 0; i < (int) (unconsMols.size()); i++)
380	<	nUnconsAtoms_local += unconsMols[i]->getNAtoms();
381	<
382	<	#ifndef IS_MPI
383	<	totNumOfUnconsAtoms = nUnconsAtoms_local;
384	<	#else
385	<	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM,
386	<	MPI_COMM_WORLD);
387	<	#endif
388	<
389	<	forcePolicy->update();
390	<	}
391	<
392	<	template<typename T> ZConstraint<T>::~ZConstraint(){
393	<
394	<	if (fzOut){
395	<	delete fzOut;
396	<	}
397	<
398	<	if (forcePolicy){
399	<	delete forcePolicy;
400	<	}
401	<	}
402	<
403	<
404	<	/**
405	<	*
406	<	*/
407	<
408	<	#ifdef IS_MPI
409	<	template<typename T> void ZConstraint<T>::update(){
410	<	double COM[3];
411	<	int index;
412	<
413	<	zconsMols.clear();
414	<	massOfZConsMols.clear();
415	<	zPos.clear();
416	<	kz.clear();
417	<	cantPos.clear();
418	<	cantVel.clear();
419	<
420	<	unconsMols.clear();
421	<	massOfUnconsMols.clear();
422	<
423	<
424	<	//creat zconsMol and unconsMol lists
425	<	for (int i = 0; i < nMols; i++){
426	<	index = isZConstraintMol(&molecules[i]);
427	<
428	<	if (index > -1){
429	<	zconsMols.push_back(&molecules[i]);
430	<	zPos.push_back((*parameters)[index].zPos);
431	<	kz.push_back((*parameters)[index].kRatio * zForceConst);
432	<	massOfZConsMols.push_back(molecules[i].getTotalMass());
433	<
434	<	if(usingSMD)
435	<	cantVel.push_back((*parameters)[index].cantVel);
436	<
437	<	}
438	<	else{
439	<	unconsMols.push_back(&molecules[i]);
440	<	massOfUnconsMols.push_back(molecules[i].getTotalMass());
441	<	}
442	<	}
443	<
444	<	fz.resize(zconsMols.size());
445	<	curZPos.resize(zconsMols.size());
446	<	indexOfZConsMols.resize(zconsMols.size());
447	<
448	<	for (size_t i = 0; i < zconsMols.size(); i++){
449	<	indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
450	<	}
451	<
452	<	//determine the states of z-constraint molecules
453	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
454	<
455	<	zconsMols[i]->getCOM(COM);
456	<
457	<	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
458	<	states.push_back(zcsFixed);
459	<
460	<	if (hasZConsGap)
461	<	endFixTime.push_back(info->getTime() + zconsFixTime);
462	<	}
463	<	else{
464	<	states.push_back(zcsMoving);
465	<
466	<	if (hasZConsGap)
467	<	endFixTime.push_back(INFINITE_TIME);
468	<	}
469	<
470	<	if(usingSMD)
471	<	cantPos.push_back(COM[whichDirection]);
472	<	}
473	<	//
474	<	forcePolicy->update();
475	<	}
476	<
477	<	#endif
478	<
479	<	/**
480	<	*  Function Name: isZConstraintMol
481	<	*  Parameter
482	<	*    Molecule* mol
483	<	*  Return value:
484	<	*    -1, if the molecule is not z-constraint molecule,
485	<	*    other non-negative values, its index in indexOfAllZConsMols vector
486	<	*/
487	<
488	<	template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
489	<	int index;
490	<	int low;
491	<	int high;
492	<	int mid;
493	<
494	<	index = mol->getGlobalIndex();
495	<
496	<	low = 0;
497	<	high = parameters->size() - 1;
498	<
499	<	//Binary Search (we have sorted the array)
500	<	while (low <= high){
501	<	mid = (low + high) / 2;
502	<	if ((*parameters)[mid].zconsIndex == index)
503	<	return mid;
504	<	else if ((*parameters)[mid].zconsIndex > index)
505	<	high = mid - 1;
506	<	else
507	<	low = mid + 1;
508	<	}
509	<
510	<	return -1;
511	<	}
512	<
513	<	template<typename T> void ZConstraint<T>::integrate(){
514	<	// creat zconsWriter
515	<	fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
516	<
517	<	if (!fzOut){
518	<	sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
519	<	painCave.isFatal = 1;
520	<	simError();
521	<	}
522	<
523	<	//zero out the velocities of center of mass of unconstrained molecules
524	<	//and the velocities of center of mass of every single z-constrained molecueles
525	<	zeroOutVel();
526	<
527	<	curZconsTime = zconsTime + info->getTime();
528	<
529	<	T::integrate();
530	<	}
531	<
532	<
533	<	/**
534	<	*
535	<	*
536	<	*
537	<	*
538	<	*/
539	<	template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
540	<	double zsys;
541	<	double COM[3];
542	<	double force[3];
543	<	double zSysCOMVel;
544	<
545	<	T::calcForce(calcPot, calcStress);
546	<
547	<
548	<	if (hasZConsGap){
549	<	updateZPos();
550	<	}
551	<
552	<	if (checkZConsState()){
553	<	zeroOutVel();
554	<	forcePolicy->update();
555	<	}
556	<
557	<	zsys = calcZSys();
558	<	zSysCOMVel = calcSysCOMVel();
559	<	#ifdef IS_MPI
560	<	if (worldRank == 0){
561	<	#endif
562	<	//cout << "---------------------------------------------------------------------" <<endl;
563	<	//cout << "current time: " << info->getTime() << endl;
564	<	//cout << "center of mass at z: " << zsys << endl;
565	<	//cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
566	<
567	<	#ifdef IS_MPI
568	<	}
569	<	#endif
570	<
571	<	//do zconstraint force;
572	<	if (haveFixedZMols()){
573	<	this->doZconstraintForce();
574	<	}
575	<
576	<	//use external force to move the molecules to the specified positions
577	<	if (haveMovingZMols()){
578	<	if (usingSMD)
579	<	this->doHarmonic(cantPos);
580	<	else
581	<	this->doHarmonic(zPos);
582	<	}
583	<
584	<	//write out forces and current positions of z-constraint molecules
585	<	if (info->getTime() >= curZconsTime){
586	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
587	<	zconsMols[i]->getCOM(COM);
588	<	curZPos[i] = COM[whichDirection];
589	<
590	<	//if the z-constraint molecule is still moving, just record its force
591	<	if (states[i] == zcsMoving){
592	<	fz[i] = 0;
593	<	Atom** movingZAtoms;
594	<	movingZAtoms = zconsMols[i]->getMyAtoms();
595	<	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
596	<	movingZAtoms[j]->getFrc(force);
597	<	fz[i] += force[whichDirection];
598	<	}
599	<	}
600	<	}
601	<	fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
602	<	&curZPos[0], &zPos[0]);
603	<	curZconsTime += zconsTime;
604	<	}
605	<
606	<	zSysCOMVel = calcSysCOMVel();
607	<	#ifdef IS_MPI
608	<	if (worldRank == 0){
609	<	#endif
610	<	//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
611	<	#ifdef IS_MPI
612	<	}
613	<	#endif
614	<	}
615	<
616	<
617	<	/**
618	<	*
619	<	*/
620	<
621	<	template<typename T> double ZConstraint<T>::calcZSys(){
622	<	//calculate reference z coordinate for z-constraint molecules
623	<	double totalMass_local;
624	<	double totalMass;
625	<	double totalMZ_local;
626	<	double totalMZ;
627	<	double massOfCurMol;
628	<	double COM[3];
629	<
630	<	totalMass_local = 0;
631	<	totalMZ_local = 0;
632	<
633	<	for (int i = 0; i < nMols; i++){
634	<	massOfCurMol = molecules[i].getTotalMass();
635	<	molecules[i].getCOM(COM);
636	<
637	<	totalMass_local += massOfCurMol;
638	<	totalMZ_local += massOfCurMol * COM[whichDirection];
639	<	}
640	<
641	<
642	<	#ifdef IS_MPI
643	<	MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM,
644	<	MPI_COMM_WORLD);
645	<	MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
646	<	#else
647	<	totalMass = totalMass_local;
648	<	totalMZ = totalMZ_local;
649	<	#endif
650	<
651	<	double zsys;
652	<	zsys = totalMZ / totalMass;
653	<
654	<	return zsys;
655	<	}
656	<
657	<	/**
658	<	*
659	<	*/
660	<	template<typename T> void ZConstraint<T>::thermalize(void){
661	<	T::thermalize();
662	<	zeroOutVel();
663	<	}
664	<
665	<	/**
666	<	*
667	<	*/
668	<
669	<	template<typename T> void ZConstraint<T>::zeroOutVel(){
670	<	Atom** fixedZAtoms;
671	<	double COMvel[3];
672	<	double vel[3];
673	<	double zSysCOMVel;
674	<
675	<	//zero out the velocities of center of mass of fixed z-constrained molecules
676	<
677	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
678	<	if (states[i] == zcsFixed){
679	<	zconsMols[i]->getCOMvel(COMvel);
680	<	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
681	<
682	<	fixedZAtoms = zconsMols[i]->getMyAtoms();
683	<
684	<	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
685	<	fixedZAtoms[j]->getVel(vel);
686	<	vel[whichDirection] -= COMvel[whichDirection];
687	<	fixedZAtoms[j]->setVel(vel);
688	<	}
689	<
690	<	zconsMols[i]->getCOMvel(COMvel);
691	<	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
692	<	}
693	<	}
694	<
695	<	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
696	<
697	<	zSysCOMVel = calcSysCOMVel();
698	<	#ifdef IS_MPI
699	<	if (worldRank == 0){
700	<	#endif
701	<	//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
702	<	#ifdef IS_MPI
703	<	}
704	<	#endif
705	<
706	<	// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
707	<	double MVzOfMovingMols_local;
708	<	double MVzOfMovingMols;
709	<	double totalMassOfMovingZMols_local;
710	<	double totalMassOfMovingZMols;
711	<
712	<	MVzOfMovingMols_local = 0;
713	<	totalMassOfMovingZMols_local = 0;
714	<
715	<	for (int i = 0; i < (int) (unconsMols.size()); i++){
716	<	unconsMols[i]->getCOMvel(COMvel);
717	<	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
718	<	}
719	<
720	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
721	<	if (states[i] == zcsMoving){
722	<	zconsMols[i]->getCOMvel(COMvel);
723	<	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
724	<	totalMassOfMovingZMols_local += massOfZConsMols[i];
725	<	}
726	<	}
727	<
728	<	#ifndef IS_MPI
729	<	MVzOfMovingMols = MVzOfMovingMols_local;
730	<	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
731	<	#else
732	<	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
733	<	MPI_SUM, MPI_COMM_WORLD);
734	<	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
735	<	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
736	<	#endif
737	<
738	<	double vzOfMovingMols;
739	<	vzOfMovingMols = MVzOfMovingMols /
740	<	(totalMassOfUncons + totalMassOfMovingZMols);
741	<
742	<	//modify the velocites of unconstrained molecules
743	<	Atom** unconsAtoms;
744	<	for (int i = 0; i < (int) (unconsMols.size()); i++){
745	<	unconsAtoms = unconsMols[i]->getMyAtoms();
746	<	for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
747	<	unconsAtoms[j]->getVel(vel);
748	<	vel[whichDirection] -= vzOfMovingMols;
749	<	unconsAtoms[j]->setVel(vel);
750	<	}
751	<	}
752	<
753	<	//modify the velocities of moving z-constrained molecuels
754	<	Atom** movingZAtoms;
755	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
756	<	if (states[i] == zcsMoving){
757	<	movingZAtoms = zconsMols[i]->getMyAtoms();
758	<	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
759	<	movingZAtoms[j]->getVel(vel);
760	<	vel[whichDirection] -= vzOfMovingMols;
761	<	movingZAtoms[j]->setVel(vel);
762	<	}
763	<	}
764	<	}
765	<
766	<
767	<	zSysCOMVel = calcSysCOMVel();
768	<	#ifdef IS_MPI
769	<	if (worldRank == 0){
770	<	#endif
771	<	//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
772	<	#ifdef IS_MPI
773	<	}
774	<	#endif
775	<	}
776	<
777	<	/**
778	<	*
779	<	*/
780	<
781	<	template<typename T> void ZConstraint<T>::doZconstraintForce(){
782	<	Atom** zconsAtoms;
783	<	double totalFZ;
784	<	double totalFZ_local;
785	<	double COM[3];
786	<	double force[3];
787	<
788	<	//constrain the molecules which do not reach the specified positions
789	<
790	<	//Zero Out the force of z-contrained molecules
791	<	totalFZ_local = 0;
792	<
793	<	//calculate the total z-contrained force of fixed z-contrained molecules
794	<
795	<	//cout << "before zero out z-constraint force on fixed z-constraint molecuels "
796	<	//       << "total force is " << calcTotalForce() << endl;
797	<
798	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
799	<	if (states[i] == zcsFixed){
800	<	zconsMols[i]->getCOM(COM);
801	<	zconsAtoms = zconsMols[i]->getMyAtoms();
802	<
803	<	fz[i] = 0;
804	<	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
805	<	zconsAtoms[j]->getFrc(force);
806	<	fz[i] += force[whichDirection];
807	<	}
808	<	totalFZ_local += fz[i];
809	<
810	<	//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
811	<	//      <<"\tcurrent zpos: " << COM[whichDirection]
812	<	//      << "\tcurrent fz: " <<fz[i] << endl;
813	<	}
814	<	}
815	<
816	<	//calculate total z-constraint force
817	<	#ifdef IS_MPI
818	<	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
819	<	#else
820	<	totalFZ = totalFZ_local;
821	<	#endif
822	<
823	<
824	<	// apply negative to fixed z-constrained molecues;
825	<	force[0] = 0;
826	<	force[1] = 0;
827	<	force[2] = 0;
828	<
829	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
830	<	if (states[i] == zcsFixed){
831	<	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
832	<	zconsAtoms = zconsMols[i]->getMyAtoms();
833	<
834	<	for (int j = 0; j < nAtomOfCurZConsMol; j++){
835	<	//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
836	<	force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i],
837	<	zconsAtoms[j],
838	<	fz[i]);
839	<	zconsAtoms[j]->addFrc(force);
840	<	}
841	<	}
842	<	}
843	<
844	<	//cout << "after zero out z-constraint force on fixed z-constraint molecuels "
845	<	//      << "total force is " << calcTotalForce() << endl;
846	<
847	<
848	<	force[0] = 0;
849	<	force[1] = 0;
850	<	force[2] = 0;
851	<
852	<	//modify the forces of unconstrained molecules
853	<	for (int i = 0; i < (int) (unconsMols.size()); i++){
854	<	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
855	<
856	<	for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
857	<	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
858	<	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],
859	<	totalFZ);
860	<	unconsAtoms[j]->addFrc(force);
861	<	}
862	<	}
863	<
864	<	//modify the forces of moving z-constrained molecules
865	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
866	<	if (states[i] == zcsMoving){
867	<	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
868	<
869	<	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
870	<	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
871	<	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],
872	<	totalFZ);
873	<	movingZAtoms[j]->addFrc(force);
874	<	}
875	<	}
876	<	}
877	<	//  cout << "after substracting z-constraint force from moving molecuels "
878	<	//        << "total force is " << calcTotalForce()  << endl;
879	<	}
880	<
881	<	/**
882	<	*
883	<	*
884	<	*/
885	<
886	<	template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
887	<	double force[3];
888	<	double harmonicU;
889	<	double harmonicF;
890	<	double COM[3];
891	<	double diff;
892	<	double totalFZ_local;
893	<	double totalFZ;
894	<
895	<	force[0] = 0;
896	<	force[1] = 0;
897	<	force[2] = 0;
898	<
899	<	totalFZ_local = 0;
900	<
901	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
902	<	if (states[i] == zcsMoving){
903	<	zconsMols[i]->getCOM(COM);
904	<	//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
905	<	//     << "\tcurrent zpos: " << COM[whichDirection] << endl;
906	<
907	<	diff = COM[whichDirection] - resPos[i];
908	<
909	<	harmonicU = 0.5 * kz[i] * diff * diff;
910	<	info->lrPot += harmonicU;
911	<
912	<	harmonicF = -kz[i] * diff;
913	<	totalFZ_local += harmonicF;
914	<
915	<	//adjust force
916	<
917	<	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
918	<
919	<	for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
920	<	//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
921	<	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i],
922	<	movingZAtoms[j],
923	<	harmonicF);
924	<	movingZAtoms[j]->addFrc(force);
925	<	}
926	<	}
927	<	}
928	<
929	<	#ifndef IS_MPI
930	<	totalFZ = totalFZ_local;
931	<	#else
932	<	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
933	<	#endif
934	<
935	<	//cout << "before substracting harmonic force from moving molecuels "
936	<	//      << "total force is " << calcTotalForce()  << endl;
937	<
938	<	force[0] = 0;
939	<	force[1] = 0;
940	<	force[2] = 0;
941	<
942	<	//modify the forces of unconstrained molecules
943	<	for (int i = 0; i < (int) (unconsMols.size()); i++){
944	<	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
945	<
946	<	for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
947	<	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
948	<	force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j],
949	<	totalFZ);
950	<	unconsAtoms[j]->addFrc(force);
951	<	}
952	<	}
953	<
954	<	//cout << "after substracting harmonic force from moving molecuels "
955	<	//      << "total force is " << calcTotalForce()  << endl;
956	<	}
957	<
958	<	/**
959	<	*
960	<	*/
961	<
962	<	template<typename T> bool ZConstraint<T>::checkZConsState(){
963	<	double COM[3];
964	<	double diff;
965	<
966	<	int changed_local;
967	<	int changed;
968	<
969	<	changed_local = 0;
970	<
971	<	for (int i = 0; i < (int) (zconsMols.size()); i++){
972	<	zconsMols[i]->getCOM(COM);
973	<	diff = fabs(COM[whichDirection] - zPos[i]);
974	<	if (diff <= zconsTol && states[i] == zcsMoving){
975	<	states[i] = zcsFixed;
976	<	changed_local = 1;
977	<
978	<	if (hasZConsGap)
979	<	endFixTime[i] = info->getTime() + zconsFixTime;
980	<	}
981	<	else if (diff > zconsTol && states[i] == zcsFixed){
982	<	states[i] = zcsMoving;
983	<	changed_local = 1;
984	<
985	<	if (hasZConsGap)
986	<	endFixTime[i] = INFINITE_TIME;
987	<	}
988	<	}
989	<
990	<	#ifndef IS_MPI
991	<	changed = changed_local;
992	<	#else
993	<	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
994	<	#endif
995	<
996	<	return (changed > 0);
997	<	}
998	<
999	<	template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1000	<	int havingFixed_local;
1001	<	int havingFixed;
1002	<
1003	<	havingFixed_local = 0;
1004	<
1005	<	for (int i = 0; i < (int) (zconsMols.size()); i++)
1006	<	if (states[i] == zcsFixed){
1007	<	havingFixed_local = 1;
1008	<	break;
1009	<	}
1010	<
1011	<	#ifndef IS_MPI
1012	<	havingFixed = havingFixed_local;
1013	<	#else
1014	<	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
1015	<	MPI_COMM_WORLD);
1016	<	#endif
1017	<
1018	<	return (havingFixed > 0);
1019	<	}
1020	<
1021	<
1022	<	/**
1023	<	*
1024	<	*/
1025	<	template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1026	<	int havingMoving_local;
1027	<	int havingMoving;
1028	<
1029	<	havingMoving_local = 0;
1030	<
1031	<	for (int i = 0; i < (int) (zconsMols.size()); i++)
1032	<	if (states[i] == zcsMoving){
1033	<	havingMoving_local = 1;
1034	<	break;
1035	<	}
1036	<
1037	<	#ifndef IS_MPI
1038	<	havingMoving = havingMoving_local;
1039	<	#else
1040	<	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
1041	<	MPI_COMM_WORLD);
1042	<	#endif
1043	<
1044	<	return (havingMoving > 0);
1045	<	}
1046	<
1047	<	/**
1048	<	*
1049	<	*/
1050	<
1051	<	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel(){
1052	<	double MVzOfMovingMols_local;
1053	<	double MVzOfMovingMols;
1054	<	double totalMassOfMovingZMols_local;
1055	<	double totalMassOfMovingZMols;
1056	<	double COMvel[3];
1057	<
1058	<	MVzOfMovingMols_local = 0;
1059	<	totalMassOfMovingZMols_local = 0;
1060	<
1061	<	for (int i = 0; i < unconsMols.size(); i++){
1062	<	unconsMols[i]->getCOMvel(COMvel);
1063	<	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1064	<	}
1065	<
1066	<	for (int i = 0; i < zconsMols.size(); i++){
1067	<	if (states[i] == zcsMoving){
1068	<	zconsMols[i]->getCOMvel(COMvel);
1069	<	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1070	<	totalMassOfMovingZMols_local += massOfZConsMols[i];
1071	<	}
1072	<	}
1073	<
1074	<	#ifndef IS_MPI
1075	<	MVzOfMovingMols = MVzOfMovingMols_local;
1076	<	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1077	<	#else
1078	<	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,
1079	<	MPI_SUM, MPI_COMM_WORLD);
1080	<	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1,
1081	<	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1082	<	#endif
1083	<
1084	<	double vzOfMovingMols;
1085	<	vzOfMovingMols = MVzOfMovingMols /
1086	<	(totalMassOfUncons + totalMassOfMovingZMols);
1087	<
1088	<	return vzOfMovingMols;
1089	<	}
1090	<
1091	<	/**
1092	<	*
1093	<	*/
1094	<
1095	<	template<typename T> double ZConstraint<T>::calcSysCOMVel(){
1096	<	double COMvel[3];
1097	<	double tempMVz_local;
1098	<	double tempMVz;
1099	<	double massOfZCons_local;
1100	<	double massOfZCons;
1101	<
1102	<
1103	<	tempMVz_local = 0;
1104	<
1105	<	for (int i = 0 ; i < nMols; i++){
1106	<	molecules[i].getCOMvel(COMvel);
1107	<	tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection];
1108	<	}
1109	<
1110	<	massOfZCons_local = 0;
1111	<
1112	<	for (int i = 0; i < (int) (massOfZConsMols.size()); i++){
1113	<	massOfZCons_local += massOfZConsMols[i];
1114	<	}
1115	<	#ifndef IS_MPI
1116	<	massOfZCons = massOfZCons_local;
1117	<	tempMVz = tempMVz_local;
1118	<	#else
1119	<	MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM,
1120	<	MPI_COMM_WORLD);
1121	<	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
1122	<	#endif
1123	<
1124	<	return tempMVz / (totalMassOfUncons + massOfZCons);
1125	<	}
1126	<
1127	<	/**
1128	<	*
1129	<	*/
1130	<
1131	<	template<typename T> double ZConstraint<T>::calcTotalForce(){
1132	<	double force[3];
1133	<	double totalForce_local;
1134	<	double totalForce;
1135	<
1136	<	totalForce_local = 0;
1137	<
1138	<	for (int i = 0; i < nAtoms; i++){
1139	<	atoms[i]->getFrc(force);
1140	<	totalForce_local += force[whichDirection];
1141	<	}
1142	<
1143	<	#ifndef IS_MPI
1144	<	totalForce = totalForce_local;
1145	<	#else
1146	<	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM,
1147	<	MPI_COMM_WORLD);
1148	<	#endif
1149	<
1150	<	return totalForce;
1151	<	}
1152	<
1153	<	/**
1154	<	*
1155	<	*/
1156	<
1157	<	template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1158	<	//calculate the number of atoms of moving z-constrained molecules
1159	<	int nMovingZAtoms_local;
1160	<	int nMovingZAtoms;
1161	<
1162	<	nMovingZAtoms_local = 0;
1163	<	for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1164	<	if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1165	<	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1166	<	}
1167	<
1168	<	#ifdef IS_MPI
1169	<	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM,
1170	<	MPI_COMM_WORLD);
1171	<	#else
1172	<	nMovingZAtoms = nMovingZAtoms_local;
1173	<	#endif
1174	<	totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1175	<	}
1176	<
1177	<	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol,
1178	<	Atom* atom,
1179	<	double totalForce){
1180	<	return totalForce / mol->getNAtoms();
1181	<	}
1182	<
1183	<	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom,
1184	<	double totalForce){
1185	<	return totalForce / totNumOfMovingAtoms;
1186	<	}
1187	<
1188	<	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol,
1189	<	Atom* atom,
1190	<	double totalForce){
1191	<	return totalForce / mol->getNAtoms();
1192	<	}
1193	<
1194	<	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom,
1195	<	double totalForce){
1196	<	return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1197	<	}
1198	<
1199	<	/**
1200	<	*
1201	<	*/
1202	<
1203	<	template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1204	<	//calculate the number of atoms of moving z-constrained molecules
1205	<	double massOfMovingZAtoms_local;
1206	<	double massOfMovingZAtoms;
1207	<
1208	<	massOfMovingZAtoms_local = 0;
1209	<	for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++)
1210	<	if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){
1211	<	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1212	<	}
1213	<
1214	<	#ifdef IS_MPI
1215	<	MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,
1216	<	MPI_SUM, MPI_COMM_WORLD);
1217	<	#else
1218	<	massOfMovingZAtoms = massOfMovingZAtoms_local;
1219	<	#endif
1220	<	totMassOfMovingAtoms = massOfMovingZAtoms +
1221	<	zconsIntegrator->totalMassOfUncons;
1222	<	}
1223	<
1224	<	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol,
1225	<	Atom* atom,
1226	<	double totalForce){
1227	<	return totalForce * atom->getMass() / mol->getTotalMass();
1228	<	}
1229	<
1230	<	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols(Atom* atom,
1231	<	double totalForce){
1232	<	return totalForce * atom->getMass() / totMassOfMovingAtoms;
1233	<	}
1234	<
1235	<	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol,
1236	<	Atom* atom,
1237	<	double totalForce){
1238	<	return totalForce * atom->getMass() / mol->getTotalMass();
1239	<	}
1240	<
1241	<	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom,
1242	<	double totalForce){
1243	<	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1244	<	}
1245	<
1246	<	template<typename T> void ZConstraint<T>::updateZPos(){
1247	<	double curTime;
1248	<	double COM[3];
1249	<
1250	<	curTime = info->getTime();
1251	<
1252	<	for (size_t i = 0; i < zconsMols.size(); i++){
1253	<
1254	<	if (states[i] == zcsFixed && curTime >= endFixTime[i]){
1255	<	zPos[i] += zconsGap;
1256	<
1257	<	if (usingSMD){
1258	<	zconsMols[i]->getCOM(COM);
1259	<	cantPos[i] = COM[whichDirection];
1260	<	}
1261	<
1262	<	}
1263	<
1264	<	}
1265	<
1266	<	}
1267	<
1268	<	template<typename T> void ZConstraint<T>::updateCantPos(){
1269	<	double curTime;
1270	<	double dt;
1271	<
1272	<	curTime = info->getTime();
1273	<	dt = info->dt;
1274	<
1275	<	for (size_t i = 0; i < zconsMols.size(); i++){
1276	<	if (states[i] == zcsMoving){
1277	<	cantPos[i] += cantVel[i] * dt;
1278	<	}
1279	<	}
1280	<
1281	<	}
1	>	#include "Integrator.hpp"
2	>	#include "simError.h"
3	>	#include <math.h>
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
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	>	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");
56	>	painCave.isFatal = 0;
57	>	simError();
58	>
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);
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	>	}
74	>	}
75	>	}
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();
87	>	}
88	>	else{
89	>	sampleTime = dynamic_cast<DoubleData*>(data);
90	>
91	>	if (!sampleTime){
92	>	sprintf(painCave.errMsg,
93	>	"ZConstraint error: Can not get property from SimInfo\n");
94	>	painCave.isFatal = 1;
95	>	simError();
96	>	}
97	>	else{
98	>	this->zconsTime = sampleTime->getData();
99	>	}
100	>	}
101	>
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();
110	>	}
111	>	else{
112	>	filename = dynamic_cast<StringData*>(data);
113	>
114	>	if (!filename){
115	>	sprintf(painCave.errMsg,
116	>	"ZConstraint error: Can not get property from SimInfo\n");
117	>	painCave.isFatal = 1;
118	>	simError();
119	>	}
120	>	else{
121	>	this->zconsOutput = filename->getData();
122	>	}
123	>	}
124	>
125	>	//retrieve tolerance for z-constraint molecuels
126	>	data = info->getProperty(ZCONSTOL_ID);
127	>
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	>	if (!tolerance){
137	>	sprintf(painCave.errMsg,
138	>	"ZConstraint error: Can not get property from SimInfo\n");
139	>	painCave.isFatal = 1;
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	>	}
164	>
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();
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	>	(*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;
311	>
312	>	for (int i = 0; i < nMols; i++){
313	>	searchResult = isZConstraintMol(&molecules[i]);
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);
320	>	kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
321	>
322	>	if(usingSMD)
323	>	cantVel.push_back((*parameters)[searchResult].cantVel);
324	>
325	>	}
326	>	else{
327	>	unconsMols.push_back(&molecules[i]);
328	>	massOfUnconsMols.push_back(molecules[i].getTotalMass());
329	>	}
330	>	}
331	>
332	>	fz.resize(zconsMols.size());
333	>	curZPos.resize(zconsMols.size());
334	>	indexOfZConsMols.resize(zconsMols.size());
335	>
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	>	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 (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 (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(){
412	>	double COM[3];
413	>	int index;
414	>
415	>	zconsMols.clear();
416	>	massOfZConsMols.clear();
417	>	zPos.clear();
418	>	kz.clear();
419	>	cantPos.clear();
420	>	cantVel.clear();
421	>
422	>	unconsMols.clear();
423	>	massOfUnconsMols.clear();
424	>
425	>
426	>	//creat zconsMol and unconsMol lists
427	>	for (int i = 0; i < nMols; i++){
428	>	index = isZConstraintMol(&molecules[i]);
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	>	if(usingSMD)
437	>	cantVel.push_back((*parameters)[index].cantVel);
438	>
439	>	}
440	>	else{
441	>	unconsMols.push_back(&molecules[i]);
442	>	massOfUnconsMols.push_back(molecules[i].getTotalMass());
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	>	//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 (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
460	>	states.push_back(zcsFixed);
461	>
462	>	if (hasZConsGap)
463	>	endFixTime.push_back(info->getTime() + zconsFixTime);
464	>	}
465	>	else{
466	>	states.push_back(zcsMoving);
467	>
468	>	if (hasZConsGap)
469	>	endFixTime.push_back(INFINITE_TIME);
470	>	}
471	>
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	>	/**
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){
495	>	int index;
496	>	int low;
497	>	int high;
498	>	int mid;
499	>
500	>	index = mol->getGlobalIndex();
501	>
502	>	low = 0;
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 ((*parameters)[mid].zconsIndex == index)
509	>	return mid;
510	>	else if ((*parameters)[mid].zconsIndex > index)
511	>	high = mid - 1;
512	>	else
513	>	low = mid + 1;
514	>	}
515	>
516	>	return -1;
517	>	}
518	>
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	>
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> 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	>	zsys = totalMZ / totalMass;
659	>
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	>	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	>	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 /
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	>
772	>
773	>	zSysCOMVel = calcSysCOMVel();
774	>	#ifdef IS_MPI
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	>	totalFZ = totalFZ_local;
827	>	#endif
828	>
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	>	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