src/restraints/Restraints.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
// Thermodynamic integration is not multiprocessor fristd::endly right now


#include <iostream>

#include <stdlib.h>

#include <cstdio>

#include <fstream>

#include <iomanip>

#include <string>

#include <cstring>

#include <math.h>


#include "restraints/Restraints.hpp"

#include "brains/SimInfo.hpp"

#include "utils/simError.h"


#define PI 3.14159265359

#define TWO_PI 6.28318530718


Restraints::Restraints(double lambdaVal, double lambdaExp){

  lambdaValue = lambdaVal;

  lambdaK = lambdaExp;

   std::vector<double> resConsts;

  const char *jolt = " \t\n;,";


#ifdef IS_MPI

  if(worldRank == 0 ){

#endif // is_mpi


    strcpy(springName, "HarmSpringConsts.txt");

    
    ifstream springs(springName);

    
    if (!springs) { 

      sprintf(painCave.errMsg,

              "Unable to open HarmSpringConsts.txt for reading, so the\n"

              "\tdefault spring constants will be loaded. If you want\n"

              "\tto specify spring constants, include a three line\n"

              "\tHarmSpringConsts.txt file in the execution directory.\n");

      painCave.severity = OOPSE_WARNING;

      painCave.isFatal = 0;

      simError();   

      
      // load default spring constants

      kDist  = 6;  // spring constant in units of kcal/(mol*ang^2)

      kTheta = 7.5;   // in units of kcal/mol

      kOmega = 13.5;   // in units of kcal/mol

    } else  {

      
      springs.getline(inLine,999,'\n');

      // the file is blank!

      if (springs.eof()){

      sprintf(painCave.errMsg,

              "HarmSpringConsts.txt file is not valid.\n"

              "\tThe file should contain four rows, the last three containing\n"

              "\ta label and the spring constant value. They should be listed\n"

              "\tin the following order: kDist (positional restrant), kTheta\n"

              "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"

              "\trestraint: rotation about the z-axis).\n");

      painCave.severity = OOPSE_ERROR;

      painCave.isFatal = 1;

      simError();   

      }

      // read in spring constants and check to make sure it is a valid file

      springs.getline(inLine,999,'\n');

      while (!springs.eof()){

        if (NULL != inLine){

          token = strtok(inLine,jolt);

          token = strtok(NULL,jolt);

          if (NULL != token){

            strcpy(inValue,token);

            resConsts.push_back(atof(inValue));

          }

        }

        springs.getline(inLine,999,'\n');

      }

      if (resConsts.size() == 3){

        kDist = resConsts[0];

        kTheta = resConsts[1];

        kOmega = resConsts[2];

      }

      else {

        sprintf(painCave.errMsg,

                "HarmSpringConsts.txt file is not valid.\n"

                "\tThe file should contain four rows, the last three containing\n"

                "\ta label and the spring constant value. They should be listed\n"

                "\tin the following order: kDist (positional restrant), kTheta\n"

                "\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n"

                "\trestraint: rotation about the z-axis).\n");

        painCave.severity = OOPSE_ERROR;

        painCave.isFatal = 1;

        simError();  

      }

    }

#ifdef IS_MPI

  }

  
  MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); 

  
  sprintf( checkPointMsg,

           "Sucessfully opened and read spring file.\n");

  MPIcheckPoint();


#endif // is_mpi

  
  sprintf(painCave.errMsg,

          "The spring constants for thermodynamic integration are:\n"

          "\tkDist = %lf\n"

          "\tkTheta = %lf\n"

          "\tkOmega = %lf\n", kDist, kTheta, kOmega);

  painCave.severity = OOPSE_INFO;

  painCave.isFatal = 0;

  simError();   

}


Restraints::~Restraints(){

}


void Restraints::Calc_rVal(double position[3], int currentMol){

  delRx = position[0] - cofmPosX[currentMol];

  delRy = position[1] - cofmPosY[currentMol];

  delRz = position[2] - cofmPosZ[currentMol];


  return;

}


void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){

  ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol]

    + matrix[0][2]*uZ0[currentMol];

  ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol]

    + matrix[1][2]*uZ0[currentMol];

  ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol]

    + matrix[2][2]*uZ0[currentMol];


  normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z);

  ub0x = ub0x/normalize;

  ub0y = ub0y/normalize;

  ub0z = ub0z/normalize;


  // Theta is the dot product of the reference and new z-axes

  theta = acos(ub0z);


  return;

}


void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){

  double zRotator[3][3];

  double tempOmega;

  double wholeTwoPis;

  // Use the omega accumulated from the rotation propagation

  omega = zAngle;


  // translate the omega into a range between -PI and PI

  if (omega < -PI){

    tempOmega = omega / -TWO_PI;

    wholeTwoPis = floor(tempOmega);

    tempOmega = omega + TWO_PI*wholeTwoPis;

    if (tempOmega < -PI)

      omega = tempOmega + TWO_PI;

    else

      omega = tempOmega;

  }

  if (omega > PI){

    tempOmega = omega / TWO_PI;

    wholeTwoPis = floor(tempOmega);

    tempOmega = omega - TWO_PI*wholeTwoPis;

    if (tempOmega > PI)

      omega = tempOmega - TWO_PI;   

    else

      omega = tempOmega;

  }


  vb0x = sin(omega);

  vb0y = cos(omega);

  vb0z = 0.0;


  normalize = sqrt(vb0x*vb0x + vb0y*vb0y + vb0z*vb0z);

  vb0x = vb0x/normalize;

  vb0y = vb0y/normalize;

  vb0z = vb0z/normalize;


  return;

}


double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){

  double pos[3];

  double A[3][3];

  double tolerance;

  double tempPotent;

  double factor;

  double spaceTrq[3];

  double omegaPass;


  tolerance = 5.72957795131e-7;


  harmPotent = 0.0;  // zero out the global harmonic potential variable


  factor = 1 - pow(lambdaValue, lambdaK);


  for (i=0; i<vecParticles.size(); i++){

    if (vecParticles[i]->isDirectional()){

      pos = vecParticles[i]->getPos();

      vecParticles[i]->getA(A);

      Calc_rVal( pos, i );

      Calc_body_thetaVal( A, i );

      omegaPass = vecParticles[i]->getZangle();

      Calc_body_omegaVal( A, omegaPass );


      // first we calculate the derivatives

      dVdrx = -kDist*delRx;

      dVdry = -kDist*delRy;

      dVdrz = -kDist*delRz;


      // uTx... and vTx... are the body-fixed z and y unit vectors

      uTx = 0.0;

      uTy = 0.0;

      uTz = 1.0;

      vTx = 0.0;

      vTy = 1.0;

      vTz = 0.0;


      dVdux = 0;

      dVduy = 0;

      dVduz = 0;

      dVdvx = 0;

      dVdvy = 0;

      dVdvz = 0;


      if (fabs(theta) > tolerance) {

        dVdux = -(kTheta*theta/sin(theta))*ub0x;

        dVduy = -(kTheta*theta/sin(theta))*ub0y;

        dVduz = -(kTheta*theta/sin(theta))*ub0z;

      }


      if (fabs(omega) > tolerance) {

        dVdvx = -(kOmega*omega/sin(omega))*vb0x;

        dVdvy = -(kOmega*omega/sin(omega))*vb0y;

        dVdvz = -(kOmega*omega/sin(omega))*vb0z;

      }


      // next we calculate the restraint forces and torques

      restraintFrc[0] = dVdrx;

      restraintFrc[1] = dVdry;

      restraintFrc[2] = dVdrz;

      tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz);


      restraintTrq[0] = 0.0;

      restraintTrq[1] = 0.0;

      restraintTrq[2] = 0.0;


      if (fabs(omega) > tolerance) {

        restraintTrq[0] += 0.0;

        restraintTrq[1] += 0.0;

        restraintTrq[2] += vTy*dVdvx;

        tempPotent += 0.5*(kOmega*omega*omega);

      }

      if (fabs(theta) > tolerance) {

        restraintTrq[0] += (uTz*dVduy);

        restraintTrq[1] += -(uTz*dVdux);

        restraintTrq[2] += 0.0;

        tempPotent += 0.5*(kTheta*theta*theta);

      }


      for (j = 0; j < 3; j++) {

        restraintFrc[j] *= factor;

        restraintTrq[j] *= factor;

      }


      harmPotent += tempPotent;


      // now we need to convert from body-fixed torques to space-fixed torques

      spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1] 

        + A[2][0]*restraintTrq[2];

      spaceTrq[1] = A[0][1]*restraintTrq[0] + A[1][1]*restraintTrq[1] 

        + A[2][1]*restraintTrq[2];

      spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1] 

        + A[2][2]*restraintTrq[2];


      // now it's time to pass these temporary forces and torques

      // to the total forces and torques

      vecParticles[i]->addFrc(restraintFrc);

      vecParticles[i]->addTrq(spaceTrq);

    }

  }


  // and we can return the appropriately scaled potential energy

  tempPotent = harmPotent * factor;

  return tempPotent;

}


void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){

  int idealSize;

  double pos[3];

  double A[3][3];

  double RfromQ[3][3];

  double quat0, quat1, quat2, quat3;

  double dot;

   std::vector<double> tempZangs;

  const char *delimit = " \t\n;,";


  //open the idealCrystal.in file and zAngle.ang file

  strcpy(fileName, "idealCrystal.in");

  strcpy(angleName, "zAngle.ang");

  
  ifstream crystalIn(fileName);

  ifstream angleIn(angleName);


  // check to see if these files are present in the execution directory

  if (!crystalIn) { 

    sprintf(painCave.errMsg,

            "Restraints Error: Unable to open idealCrystal.in for reading.\n"

            "\tMake sure a ref. crystal file is in the working directory.\n");

    painCave.severity = OOPSE_ERROR;

    painCave.isFatal = 1;

    simError();   

  }


  // it's not fatal to lack a zAngle.ang file, it just means you're starting

  // from the ideal crystal state

  if (!angleIn) { 

    sprintf(painCave.errMsg,

            "Restraints Warning: The lack of a zAngle.ang file is mildly\n"

            "\tunsettling... This means the simulation is starting from the\n"

            "\tidealCrystal.in reference configuration, so the omega values\n"

            "\twill all be set to zero. If this is not the case, the energy\n"

            "\tcalculations will be wrong.\n");

    painCave.severity = OOPSE_WARNING;

    painCave.isFatal = 0;

    simError();   

  }


  // A rather specific reader for OOPSE .eor files...

  // Let's read in the perfect crystal file

  crystalIn.getline(inLine,999,'\n');

  // check to see if the crystal file is the same length as starting config.

  token = strtok(inLine,delimit);

  strcpy(inValue,token);

  idealSize = atoi(inValue);

  if (idealSize != vecParticles.size()) {

    sprintf(painCave.errMsg,

            "Restraints Error: Reference crystal file is not valid.\n"

            "\tMake sure the idealCrystal.in file is the same size as the\n"

            "\tstarting configuration. Using an incompatable crystal will\n"

            "\tlead to energy calculation failures.\n");

    painCave.severity = OOPSE_ERROR;

    painCave.isFatal = 1;

    simError();  

  }

  // else, the file is okay... let's continue

  crystalIn.getline(inLine,999,'\n');

  
  for (i=0; i<vecParticles.size(); i++) {

    crystalIn.getline(inLine,999,'\n');

    token = strtok(inLine,delimit);

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    cofmPosX.push_back(atof(inValue));

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    cofmPosY.push_back(atof(inValue));

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    cofmPosZ.push_back(atof(inValue));

    token = strtok(NULL,delimit);

    token = strtok(NULL,delimit);

    token = strtok(NULL,delimit);

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    quat0 = atof(inValue);

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    quat1 = atof(inValue);

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    quat2 = atof(inValue);

    token = strtok(NULL,delimit);

    strcpy(inValue,token);

    quat3 = atof(inValue);


    // now build the rotation matrix and find the unit vectors

    RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3;

    RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3);

    RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2);

    RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3);

    RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3;

    RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1);

    RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2);

    RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1);

    RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3;

    
    normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] 

                     + RfromQ[2][2]*RfromQ[2][2]);

    uX0.push_back(RfromQ[2][0]/normalize);

    uY0.push_back(RfromQ[2][1]/normalize);

    uZ0.push_back(RfromQ[2][2]/normalize);


    normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1]

                     + RfromQ[1][2]*RfromQ[1][2]);

    vX0.push_back(RfromQ[1][0]/normalize);

    vY0.push_back(RfromQ[1][1]/normalize);

    vZ0.push_back(RfromQ[1][2]/normalize);

  }

  crystalIn.close();


  // now we read in the zAngle.ang file

  if (angleIn){

    angleIn.getline(inLine,999,'\n');

    angleIn.getline(inLine,999,'\n');

    while (!angleIn.eof()) {

      token = strtok(inLine,delimit);

      strcpy(inValue,token);

      tempZangs.push_back(atof(inValue));

      angleIn.getline(inLine,999,'\n');

    }


    // test to make sure the zAngle.ang file is the proper length

    if (tempZangs.size() == vecParticles.size())

      for (i=0; i<vecParticles.size(); i++)

        vecParticles[i]->setZangle(tempZangs[i]);

    else {

      sprintf(painCave.errMsg,

              "Restraints Error: the supplied zAngle file is not valid.\n"

              "\tMake sure the zAngle.ang file matches with the initial\n"

              "\tconfiguration (i.e. they're the same length). Using the wrong\n"

              "\tzAngle file will lead to errors in the energy calculations.\n");

      painCave.severity = OOPSE_ERROR;

      painCave.isFatal = 1;

      simError();  

    }

  }

  angleIn.close();


  return;

}


void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){


  char zOutName[200];


  strcpy(zOutName,"zAngle.ang");


  ofstream angleOut(zOutName);

  angleOut << "This file contains the omega values for the .eor file\n";

  for (i=0; i<vecParticles.size(); i++) {

    angleOut << vecParticles[i]->getZangle() << "\n";

  }

  return;

}


double Restraints::getVharm(){

  return harmPotent;

}


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