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 733 by tim, Wed Aug 27 19:23:29 2003 UTC vs.
Revision 1093 by tim, Wed Mar 17 14:22:59 2004 UTC

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

Diff Legend

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