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 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/restraints/Restraints.cpp (file contents), Revision 1826 by tim, Thu Dec 2 05:04:20 2004 UTC

# Line 1 | Line 1
1 + // Thermodynamic integration is not multiprocessor friendly 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  
21 +
22 +
23   using namespace std;
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  
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 <              "In Restraints: Unable to open HarmSpringConsts.txt for reading.\n"
73 <              "\tDefault spring constants will be loaded. If you want to specify\n"
74 <              "\tspring constants, include a three line HarmSpringConsts.txt file\n"
75 <              "\tin the current directory.\n");
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 <      token = strtok(inLine,jolt);
133 <      token = strtok(NULL,jolt);
134 <      strcpy(inValue,token);
135 <      kDist = (atof(inValue));
136 <      springs.getline(inLine,999,'\n');
137 <      token = strtok(inLine,jolt);
138 <      token = strtok(NULL,jolt);
139 <      strcpy(inValue,token);
140 <      kTheta = (atof(inValue));
141 <      springs.getline(inLine,999,'\n');
142 <      token = strtok(inLine,jolt);
143 <      token = strtok(NULL,jolt);
144 <      strcpy(inValue,token);
145 <      kOmega = (atof(inValue));
146 <      springs.close();
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  
183 –      if (omega > PI || omega < -PI)
184 –        cout << "oops... " << omega << "\n";
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 < //   char *token;
601 < //   char fileName[200];
602 < //   char angleName[200];
277 < //   char inLine[1000];
278 < //   char inValue[200];
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 <            "\tMake sure a reference crystal file is in the current directory.\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 <    
295 <    return;
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 <            "\twill all be set to zero. If this is not the case, you should\n"
655 <            "\tquestion your results.\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 <  // now we can read in the zAngle.ang file
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 <    for (i=0; i<vecParticles.size(); i++) {
821 <      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 <      vecParticles[i]->setZangle(atof(inValue));
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