ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/restraints/Restraints.cpp
(Generate patch)

Comparing:
trunk/OOPSE-3.0/src/restraints/Restraints.cpp (file contents), Revision 1537 by gezelter, Wed Oct 6 21:27:21 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/restraints/Restraints.cpp (file contents), Revision 1830 by tim, Thu Dec 2 05:17:10 2004 UTC

# Line 1 | Line 1
1 < // Thermodynamic integration is not multiprocessor friendly right now
1 > // Thermodynamic integration is not multiprocessor fristd::endly right now
2  
3 +
4 +
5   #include <iostream>
6 +
7   #include <stdlib.h>
8 +
9   #include <cstdio>
10 +
11   #include <fstream>
12 +
13   #include <iomanip>
14 +
15   #include <string>
16 +
17   #include <cstring>
18 +
19   #include <math.h>
20  
12 – using namespace std;
21  
22 +
23 +
24 +
25 +
26 +
27   #include "restraints/Restraints.hpp"
28 +
29   #include "brains/SimInfo.hpp"
30 +
31   #include "utils/simError.h"
32  
33 +
34 +
35   #define PI 3.14159265359
36 +
37   #define TWO_PI 6.28318530718
38  
39 +
40 +
41   Restraints::Restraints(double lambdaVal, double lambdaExp){
42 +
43    lambdaValue = lambdaVal;
44 +
45    lambdaK = lambdaExp;
46 <  vector<double> resConsts;
46 >
47 >   std::vector<double> resConsts;
48 >
49    const char *jolt = " \t\n;,";
50  
51 +
52 +
53   #ifdef IS_MPI
54 +
55    if(worldRank == 0 ){
56 +
57   #endif // is_mpi
58  
59 +
60 +
61      strcpy(springName, "HarmSpringConsts.txt");
62 +
63      
64 +
65      ifstream springs(springName);
66 +
67      
68 +
69      if (!springs) {
70 +
71        sprintf(painCave.errMsg,
72 +
73                "Unable to open HarmSpringConsts.txt for reading, so the\n"
74 +
75                "\tdefault spring constants will be loaded. If you want\n"
76 +
77                "\tto specify spring constants, include a three line\n"
78 +
79                "\tHarmSpringConsts.txt file in the execution directory.\n");
80 +
81        painCave.severity = OOPSE_WARNING;
82 +
83        painCave.isFatal = 0;
84 +
85        simError();  
86 +
87        
88 +
89        // load default spring constants
90 +
91        kDist  = 6;  // spring constant in units of kcal/(mol*ang^2)
92 +
93        kTheta = 7.5;   // in units of kcal/mol
94 +
95        kOmega = 13.5;   // in units of kcal/mol
96 +
97      } else  {
98 +
99        
100 +
101        springs.getline(inLine,999,'\n');
102 +
103        // the file is blank!
104 +
105        if (springs.eof()){
106 +
107        sprintf(painCave.errMsg,
108 +
109                "HarmSpringConsts.txt file is not valid.\n"
110 +
111                "\tThe file should contain four rows, the last three containing\n"
112 +
113                "\ta label and the spring constant value. They should be listed\n"
114 +
115                "\tin the following order: kDist (positional restrant), kTheta\n"
116 +
117                "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
118 +
119                "\trestraint: rotation about the z-axis).\n");
120 +
121        painCave.severity = OOPSE_ERROR;
122 +
123        painCave.isFatal = 1;
124 +
125        simError();  
126 +
127        }
128 +
129        // read in spring constants and check to make sure it is a valid file
130 +
131        springs.getline(inLine,999,'\n');
132 +
133        while (!springs.eof()){
134 +
135          if (NULL != inLine){
136 +
137            token = strtok(inLine,jolt);
138 +
139            token = strtok(NULL,jolt);
140 +
141            if (NULL != token){
142 +
143              strcpy(inValue,token);
144 +
145              resConsts.push_back(atof(inValue));
146 +
147            }
148 +
149          }
150 +
151          springs.getline(inLine,999,'\n');
152 +
153        }
154 +
155        if (resConsts.size() == 3){
156 +
157          kDist = resConsts[0];
158 +
159          kTheta = resConsts[1];
160 +
161          kOmega = resConsts[2];
162 +
163        }
164 +
165        else {
166 +
167          sprintf(painCave.errMsg,
168 +
169                  "HarmSpringConsts.txt file is not valid.\n"
170 +
171                  "\tThe file should contain four rows, the last three containing\n"
172 +
173                  "\ta label and the spring constant value. They should be listed\n"
174 +
175                  "\tin the following order: kDist (positional restrant), kTheta\n"
176 +
177                  "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"
178 +
179                  "\trestraint: rotation about the z-axis).\n");
180 +
181          painCave.severity = OOPSE_ERROR;
182 +
183          painCave.isFatal = 1;
184 +
185          simError();  
186 +
187        }
188 +
189      }
190 +
191   #ifdef IS_MPI
192 +
193    }
194 +
195    
196 +
197    MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
198 +
199    MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
200 +
201    MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
202 +
203    
204 +
205    sprintf( checkPointMsg,
206 +
207             "Sucessfully opened and read spring file.\n");
208 +
209    MPIcheckPoint();
210  
211 +
212 +
213   #endif // is_mpi
214 +
215    
216 +
217    sprintf(painCave.errMsg,
218 +
219            "The spring constants for thermodynamic integration are:\n"
220 +
221            "\tkDist = %lf\n"
222 +
223            "\tkTheta = %lf\n"
224 +
225            "\tkOmega = %lf\n", kDist, kTheta, kOmega);
226 +
227    painCave.severity = OOPSE_INFO;
228 +
229    painCave.isFatal = 0;
230 +
231    simError();  
232 +
233   }
234  
235 +
236 +
237   Restraints::~Restraints(){
238 +
239   }
240  
241 +
242 +
243   void Restraints::Calc_rVal(double position[3], int currentMol){
244 +
245    delRx = position[0] - cofmPosX[currentMol];
246 +
247    delRy = position[1] - cofmPosY[currentMol];
248 +
249    delRz = position[2] - cofmPosZ[currentMol];
250  
251 +
252 +
253    return;
254 +
255   }
256  
257 +
258 +
259   void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){
260 +
261    ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol]
262 +
263      + matrix[0][2]*uZ0[currentMol];
264 +
265    ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol]
266 +
267      + matrix[1][2]*uZ0[currentMol];
268 +
269    ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol]
270 +
271      + matrix[2][2]*uZ0[currentMol];
272  
273 +
274 +
275    normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z);
276 +
277    ub0x = ub0x/normalize;
278 +
279    ub0y = ub0y/normalize;
280 +
281    ub0z = ub0z/normalize;
282  
283 +
284 +
285    // Theta is the dot product of the reference and new z-axes
286 +
287    theta = acos(ub0z);
288  
289 +
290 +
291    return;
292 +
293   }
294  
295 +
296 +
297   void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){
298 +
299    double zRotator[3][3];
300 +
301    double tempOmega;
302 +
303    double wholeTwoPis;
304 +
305    // Use the omega accumulated from the rotation propagation
306 +
307    omega = zAngle;
308  
309 +
310 +
311    // translate the omega into a range between -PI and PI
312 +
313    if (omega < -PI){
314 +
315      tempOmega = omega / -TWO_PI;
316 +
317      wholeTwoPis = floor(tempOmega);
318 +
319      tempOmega = omega + TWO_PI*wholeTwoPis;
320 +
321      if (tempOmega < -PI)
322 +
323        omega = tempOmega + TWO_PI;
324 +
325      else
326 +
327        omega = tempOmega;
328 +
329    }
330 +
331    if (omega > PI){
332 +
333      tempOmega = omega / TWO_PI;
334 +
335      wholeTwoPis = floor(tempOmega);
336 +
337      tempOmega = omega - TWO_PI*wholeTwoPis;
338 +
339      if (tempOmega > PI)
340 +
341        omega = tempOmega - TWO_PI;  
342 +
343      else
344 +
345        omega = tempOmega;
346 +
347    }
348  
349 +
350 +
351    vb0x = sin(omega);
352 +
353    vb0y = cos(omega);
354 +
355    vb0z = 0.0;
356  
357 +
358 +
359    normalize = sqrt(vb0x*vb0x + vb0y*vb0y + vb0z*vb0z);
360 +
361    vb0x = vb0x/normalize;
362 +
363    vb0y = vb0y/normalize;
364 +
365    vb0z = vb0z/normalize;
366  
367 +
368 +
369    return;
370 +
371   }
372  
373 +
374 +
375   double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){
376 +
377    double pos[3];
378 +
379    double A[3][3];
380 +
381    double tolerance;
382 +
383    double tempPotent;
384 +
385    double factor;
386 +
387    double spaceTrq[3];
388 +
389    double omegaPass;
390  
391 +
392 +
393    tolerance = 5.72957795131e-7;
394  
395 +
396 +
397    harmPotent = 0.0;  // zero out the global harmonic potential variable
398  
399 +
400 +
401    factor = 1 - pow(lambdaValue, lambdaK);
402  
403 +
404 +
405    for (i=0; i<vecParticles.size(); i++){
406 +
407      if (vecParticles[i]->isDirectional()){
408 <      vecParticles[i]->getPos(pos);
408 >
409 >      pos = vecParticles[i]->getPos();
410 >
411        vecParticles[i]->getA(A);
412 +
413        Calc_rVal( pos, i );
414 +
415        Calc_body_thetaVal( A, i );
416 +
417        omegaPass = vecParticles[i]->getZangle();
418 +
419        Calc_body_omegaVal( A, omegaPass );
420  
421 +
422 +
423        // first we calculate the derivatives
424 +
425        dVdrx = -kDist*delRx;
426 +
427        dVdry = -kDist*delRy;
428 +
429        dVdrz = -kDist*delRz;
430  
431 +
432 +
433        // uTx... and vTx... are the body-fixed z and y unit vectors
434 +
435        uTx = 0.0;
436 +
437        uTy = 0.0;
438 +
439        uTz = 1.0;
440 +
441        vTx = 0.0;
442 +
443        vTy = 1.0;
444 +
445        vTz = 0.0;
446  
447 +
448 +
449        dVdux = 0;
450 +
451        dVduy = 0;
452 +
453        dVduz = 0;
454 +
455        dVdvx = 0;
456 +
457        dVdvy = 0;
458 +
459        dVdvz = 0;
460  
461 +
462 +
463        if (fabs(theta) > tolerance) {
464 +
465          dVdux = -(kTheta*theta/sin(theta))*ub0x;
466 +
467          dVduy = -(kTheta*theta/sin(theta))*ub0y;
468 +
469          dVduz = -(kTheta*theta/sin(theta))*ub0z;
470 +
471        }
472  
473 +
474 +
475        if (fabs(omega) > tolerance) {
476 +
477          dVdvx = -(kOmega*omega/sin(omega))*vb0x;
478 +
479          dVdvy = -(kOmega*omega/sin(omega))*vb0y;
480 +
481          dVdvz = -(kOmega*omega/sin(omega))*vb0z;
482 +
483        }
484  
485 +
486 +
487        // next we calculate the restraint forces and torques
488 +
489        restraintFrc[0] = dVdrx;
490 +
491        restraintFrc[1] = dVdry;
492 +
493        restraintFrc[2] = dVdrz;
494 +
495        tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);
496  
497 +
498 +
499        restraintTrq[0] = 0.0;
500 +
501        restraintTrq[1] = 0.0;
502 +
503        restraintTrq[2] = 0.0;
504  
505 +
506 +
507        if (fabs(omega) > tolerance) {
508 +
509          restraintTrq[0] += 0.0;
510 +
511          restraintTrq[1] += 0.0;
512 +
513          restraintTrq[2] += vTy*dVdvx;
514 +
515          tempPotent += 0.5*(kOmega*omega*omega);
516 +
517        }
518 +
519        if (fabs(theta) > tolerance) {
520 +
521          restraintTrq[0] += (uTz*dVduy);
522 +
523          restraintTrq[1] += -(uTz*dVdux);
524 +
525          restraintTrq[2] += 0.0;
526 +
527          tempPotent += 0.5*(kTheta*theta*theta);
528 +
529        }
530  
531 +
532 +
533        for (j = 0; j < 3; j++) {
534 +
535          restraintFrc[j] *= factor;
536 +
537          restraintTrq[j] *= factor;
538 +
539        }
540  
541 +
542 +
543        harmPotent += tempPotent;
544  
545 +
546 +
547        // now we need to convert from body-fixed torques to space-fixed torques
548 +
549        spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1]
550 +
551          + A[2][0]*restraintTrq[2];
552 +
553        spaceTrq[1] = A[0][1]*restraintTrq[0] + A[1][1]*restraintTrq[1]
554 +
555          + A[2][1]*restraintTrq[2];
556 +
557        spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1]
558 +
559          + A[2][2]*restraintTrq[2];
560  
561 +
562 +
563        // now it's time to pass these temporary forces and torques
564 +
565        // to the total forces and torques
566 +
567        vecParticles[i]->addFrc(restraintFrc);
568 +
569        vecParticles[i]->addTrq(spaceTrq);
570 +
571      }
572 +
573    }
574  
575 +
576 +
577    // and we can return the appropriately scaled potential energy
578 +
579    tempPotent = harmPotent * factor;
580 +
581    return tempPotent;
582 +
583   }
584  
585 +
586 +
587   void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){
588 +
589    int idealSize;
590 +
591    double pos[3];
592 +
593    double A[3][3];
594 +
595    double RfromQ[3][3];
596 +
597    double quat0, quat1, quat2, quat3;
598 +
599    double dot;
600 <  vector<double> tempZangs;
600 >
601 >   std::vector<double> tempZangs;
602 >
603    const char *delimit = " \t\n;,";
604  
605 +
606 +
607    //open the idealCrystal.in file and zAngle.ang file
608 +
609    strcpy(fileName, "idealCrystal.in");
610 +
611    strcpy(angleName, "zAngle.ang");
612 +
613    
614 +
615    ifstream crystalIn(fileName);
616 +
617    ifstream angleIn(angleName);
618  
619 +
620 +
621    // check to see if these files are present in the execution directory
622 +
623    if (!crystalIn) {
624 +
625      sprintf(painCave.errMsg,
626 +
627              "Restraints Error: Unable to open idealCrystal.in for reading.\n"
628 +
629              "\tMake sure a ref. crystal file is in the working directory.\n");
630 +
631      painCave.severity = OOPSE_ERROR;
632 +
633      painCave.isFatal = 1;
634 +
635      simError();  
636 +
637    }
638  
639 +
640 +
641    // it's not fatal to lack a zAngle.ang file, it just means you're starting
642 +
643    // from the ideal crystal state
644 +
645    if (!angleIn) {
646 +
647      sprintf(painCave.errMsg,
648 +
649              "Restraints Warning: The lack of a zAngle.ang file is mildly\n"
650 +
651              "\tunsettling... This means the simulation is starting from the\n"
652 +
653              "\tidealCrystal.in reference configuration, so the omega values\n"
654 +
655              "\twill all be set to zero. If this is not the case, the energy\n"
656 +
657              "\tcalculations will be wrong.\n");
658 +
659      painCave.severity = OOPSE_WARNING;
660 +
661      painCave.isFatal = 0;
662 +
663      simError();  
664 +
665    }
666  
667 +
668 +
669    // A rather specific reader for OOPSE .eor files...
670 +
671    // Let's read in the perfect crystal file
672 +
673    crystalIn.getline(inLine,999,'\n');
674 +
675    // check to see if the crystal file is the same length as starting config.
676 +
677    token = strtok(inLine,delimit);
678 +
679    strcpy(inValue,token);
680 +
681    idealSize = atoi(inValue);
682 +
683    if (idealSize != vecParticles.size()) {
684 +
685      sprintf(painCave.errMsg,
686 +
687              "Restraints Error: Reference crystal file is not valid.\n"
688 +
689              "\tMake sure the idealCrystal.in file is the same size as the\n"
690 +
691              "\tstarting configuration. Using an incompatable crystal will\n"
692 +
693              "\tlead to energy calculation failures.\n");
694 +
695      painCave.severity = OOPSE_ERROR;
696 +
697      painCave.isFatal = 1;
698 +
699      simError();  
700 +
701    }
702 +
703    // else, the file is okay... let's continue
704 +
705    crystalIn.getline(inLine,999,'\n');
706 +
707    
708 +
709    for (i=0; i<vecParticles.size(); i++) {
710 +
711      crystalIn.getline(inLine,999,'\n');
712 +
713      token = strtok(inLine,delimit);
714 +
715      token = strtok(NULL,delimit);
716 +
717      strcpy(inValue,token);
718 +
719      cofmPosX.push_back(atof(inValue));
720 +
721      token = strtok(NULL,delimit);
722 +
723      strcpy(inValue,token);
724 +
725      cofmPosY.push_back(atof(inValue));
726 +
727      token = strtok(NULL,delimit);
728 +
729      strcpy(inValue,token);
730 +
731      cofmPosZ.push_back(atof(inValue));
732 +
733      token = strtok(NULL,delimit);
734 +
735      token = strtok(NULL,delimit);
736 +
737      token = strtok(NULL,delimit);
738 +
739      token = strtok(NULL,delimit);
740 +
741      strcpy(inValue,token);
742 +
743      quat0 = atof(inValue);
744 +
745      token = strtok(NULL,delimit);
746 +
747      strcpy(inValue,token);
748 +
749      quat1 = atof(inValue);
750 +
751      token = strtok(NULL,delimit);
752 +
753      strcpy(inValue,token);
754 +
755      quat2 = atof(inValue);
756 +
757      token = strtok(NULL,delimit);
758 +
759      strcpy(inValue,token);
760 +
761      quat3 = atof(inValue);
762  
763 +
764 +
765      // now build the rotation matrix and find the unit vectors
766 +
767      RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3;
768 +
769      RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3);
770 +
771      RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2);
772 +
773      RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3);
774 +
775      RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3;
776 +
777      RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1);
778 +
779      RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2);
780 +
781      RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1);
782 +
783      RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3;
784 +
785      
786 +
787      normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1]
788 +
789                       + RfromQ[2][2]*RfromQ[2][2]);
790 +
791      uX0.push_back(RfromQ[2][0]/normalize);
792 +
793      uY0.push_back(RfromQ[2][1]/normalize);
794 +
795      uZ0.push_back(RfromQ[2][2]/normalize);
796  
797 +
798 +
799      normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1]
800 +
801                       + RfromQ[1][2]*RfromQ[1][2]);
802 +
803      vX0.push_back(RfromQ[1][0]/normalize);
804 +
805      vY0.push_back(RfromQ[1][1]/normalize);
806 +
807      vZ0.push_back(RfromQ[1][2]/normalize);
808 +
809    }
810 +
811    crystalIn.close();
812  
813 +
814 +
815    // now we read in the zAngle.ang file
816 +
817    if (angleIn){
818 +
819      angleIn.getline(inLine,999,'\n');
820 +
821      angleIn.getline(inLine,999,'\n');
822 +
823      while (!angleIn.eof()) {
824 +
825        token = strtok(inLine,delimit);
826 +
827        strcpy(inValue,token);
828 +
829        tempZangs.push_back(atof(inValue));
830 +
831        angleIn.getline(inLine,999,'\n');
832 +
833      }
834  
835 +
836 +
837      // test to make sure the zAngle.ang file is the proper length
838 +
839      if (tempZangs.size() == vecParticles.size())
840 +
841        for (i=0; i<vecParticles.size(); i++)
842 +
843          vecParticles[i]->setZangle(tempZangs[i]);
844 +
845      else {
846 +
847        sprintf(painCave.errMsg,
848 +
849                "Restraints Error: the supplied zAngle file is not valid.\n"
850 +
851                "\tMake sure the zAngle.ang file matches with the initial\n"
852 +
853                "\tconfiguration (i.e. they're the same length). Using the wrong\n"
854 +
855                "\tzAngle file will lead to errors in the energy calculations.\n");
856 +
857        painCave.severity = OOPSE_ERROR;
858 +
859        painCave.isFatal = 1;
860 +
861        simError();  
862 +
863      }
864 +
865    }
866 +
867    angleIn.close();
868  
869 +
870 +
871    return;
872 +
873   }
874  
875 +
876 +
877   void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){
878  
879 +
880 +
881    char zOutName[200];
882  
883 +
884 +
885    strcpy(zOutName,"zAngle.ang");
886  
887 +
888 +
889    ofstream angleOut(zOutName);
890 +
891    angleOut << "This file contains the omega values for the .eor file\n";
892 +
893    for (i=0; i<vecParticles.size(); i++) {
894 +
895      angleOut << vecParticles[i]->getZangle() << "\n";
896 +
897    }
898 +
899    return;
900 +
901   }
902  
903 +
904 +
905   double Restraints::getVharm(){
906 +
907    return harmPotent;
908 +
909   }
910  
911 +
912 +

Diff Legend

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