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

Diff Legend

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