hdiff output

r32161/amber_potential.cu 2017-03-22 12:30:29.084253943 +0000 r32160/amber_potential.cu 2017-03-22 12:30:32.624300646 +0000
  6:   6: 
  7: #include "cost_function.h"  7: #include "cost_function.h"
  8: #include "potential.h"  8: #include "potential.h"
  9: #include "gpu.h"  9: #include "gpu.h"
 10:  10: 
 11: AmberPotential::AmberPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas,  11: AmberPotential::AmberPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, 
 12:                 size_t numDimensions, int nDegFreedom, int nRigidBody, int rigidMaxSite,  12:                 size_t numDimensions, int nDegFreedom, int nRigidBody, int rigidMaxSite, 
 13:                 int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody, int *rigidSingles,  13:                 int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody, int *rigidSingles, 
 14:                 double *rigidInverse, double *x, int nSecDiag, bool isAtomisticNotRigid,  14:                 double *rigidInverse, double *x, int nSecDiag, bool isAtomisticNotRigid, 
 15:                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen,  15:                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, 
 16:                 int nFreeze, bool isAaConvergence, int nAddTarget, double *ljAddRep, double *ljAddAtt) 16:                 int nFreeze, bool isAaConvergence)
 17:         : CostFunction(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite,  17:         : CostFunction(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 
 18:                         nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, x, nSecDiag,  18:                         nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, x, nSecDiag, 
 19:                         isAtomisticNotRigid, aaConvThreshold, coldFusionLim, shouldFreeze, isAtomFrozen, nFreeze,  19:                         isAtomisticNotRigid, aaConvThreshold, coldFusionLim, shouldFreeze, isAtomFrozen, nFreeze, 
 20:                         isAaConvergence, nAddTarget, ljAddRep, ljAddAtt) {} 20:                         isAaConvergence) {}
 21:  21: 
 22:  22: 
 23:  23: 
 24: void AmberPotential::computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf) 24: void AmberPotential::computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf)
 25: { 25: {
 26:         double testEnergy; 26:         double testEnergy;
 27:  27: 
 28:         // Copy coordinates to correct location in device memory.  28:         // Copy coordinates to correct location in device memory. 
 29:         gminoptim_copy_crd(d_x); 29:         gminoptim_copy_crd(d_x);
 30:         CudaCheckError(); 30:         CudaCheckError();


r32161/cost_function.cu 2017-03-22 12:30:27.984239431 +0000 r32160/cost_function.cu 2017-03-22 12:30:31.516286028 +0000
  5:  **/  5:  **/
  6:   6: 
  7: #include<iomanip>  7: #include<iomanip>
  8:   8: 
  9: #include "cost_function.h"  9: #include "cost_function.h"
 10:  10: 
 11: CostFunction::CostFunction(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom,  11: CostFunction::CostFunction(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom, 
 12:                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody,  12:                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody, 
 13:                 int *rigidSingles, double *rigidInverse, double *coords, int nSecDiag, bool isAtomisticNotRigid,  13:                 int *rigidSingles, double *rigidInverse, double *coords, int nSecDiag, bool isAtomisticNotRigid, 
 14:                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze,  14:                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze, 
 15:                 bool isAaConvergence, int nAddTarget, double *ljAddRep, double *ljAddAtt) 15:                 bool isAaConvergence)
 16:         : m_timer_potential(timer_potential) 16:         : m_timer_potential(timer_potential)
 17:         , m_numDimensions(numDimensions) 17:         , m_numDimensions(numDimensions)
 18:         , m_nDegFreedom(nDegFreedom) 18:         , m_nDegFreedom(nDegFreedom)
 19:         , m_nRigidBody(nRigidBody) 19:         , m_nRigidBody(nRigidBody)
 20:         , m_rigidMaxSite(rigidMaxSite) 20:         , m_rigidMaxSite(rigidMaxSite)
 21:         , m_nRigidSitesPerBody(nRigidSitesPerBody) 21:         , m_nRigidSitesPerBody(nRigidSitesPerBody)
 22:         , m_rigidGroups(rigidGroups) 22:         , m_rigidGroups(rigidGroups)
 23:         , m_sitesRigidBody(sitesRigidBody) 23:         , m_sitesRigidBody(sitesRigidBody)
 24:         , m_rigidSingles(rigidSingles) 24:         , m_rigidSingles(rigidSingles)
 25:         , m_rigidInverse(rigidInverse) 25:         , m_rigidInverse(rigidInverse)
 40:         , m_aaConvThreshold(aaConvThreshold) 40:         , m_aaConvThreshold(aaConvThreshold)
 41:         , m_nCalls(0) 41:         , m_nCalls(0)
 42:         , m_coldFusionLim(coldFusionLim) 42:         , m_coldFusionLim(coldFusionLim)
 43:         , m_shouldFreeze(shouldFreeze) 43:         , m_shouldFreeze(shouldFreeze)
 44:         , m_isAtomFrozen(isAtomFrozen) 44:         , m_isAtomFrozen(isAtomFrozen)
 45:         , m_nFreeze(nFreeze) 45:         , m_nFreeze(nFreeze)
 46:         , m_nLargeRigidBody(0) 46:         , m_nLargeRigidBody(0)
 47:         , m_debugPrinting(debugPrinting) 47:         , m_debugPrinting(debugPrinting)
 48:         , m_cublas(cublas) 48:         , m_cublas(cublas)
 49:         , m_isBfgsts(false) 49:         , m_isBfgsts(false)
 50:         , m_isAaConvergence(isAaConvergence) 50:           , m_isAaConvergence(isAaConvergence)
 51:         , m_nAddTarget(nAddTarget) 
 52:         , m_ljAddRep(ljAddRep) 
 53:           , m_ljAddAtt(ljAddAtt) 
 54: { 51: {
 55:         for (int i = 0; i < m_nRigidBody; ++i) { 52:         for (int i = 0; i < m_nRigidBody; ++i) {
 56:                 if (m_nRigidSitesPerBody[i] > 32) { 53:                 if (m_nRigidSitesPerBody[i] > 32) {
 57:                         m_nLargeRigidBody += 1; 54:                         m_nLargeRigidBody += 1;
 58:                 } 55:                 }
 59:         } 56:         }
 60:  57: 
 61:         m_largeRigidIndices = new int[m_nLargeRigidBody]; 58:         m_largeRigidIndices = new int[m_nLargeRigidBody];
 62:  59: 
 63:         double *zerosx; 60:         double *zerosx;
 92:                 CudaSafeCall( cudaMalloc(&m_d_coords,   numDimensions * sizeof(double)) ); 89:                 CudaSafeCall( cudaMalloc(&m_d_coords,   numDimensions * sizeof(double)) );
 93:                 CudaSafeCall( cudaMalloc(&m_d_energy, sizeof(double)) ); 90:                 CudaSafeCall( cudaMalloc(&m_d_energy, sizeof(double)) );
 94:                 CudaSafeCall( cudaMalloc(&m_d_gradient, numDimensions * sizeof(double)) ); 91:                 CudaSafeCall( cudaMalloc(&m_d_gradient, numDimensions * sizeof(double)) );
 95:         } 92:         }
 96:  93: 
 97:         CudaSafeCall( cudaMalloc(&m_d_isAtomFrozen, numDimensions/3 * sizeof(bool)) ); 94:         CudaSafeCall( cudaMalloc(&m_d_isAtomFrozen, numDimensions/3 * sizeof(bool)) );
 98:  95: 
 99:         CudaSafeCall( cudaMalloc(&m_d_nLargeRigidBody, sizeof(int)) ); 96:         CudaSafeCall( cudaMalloc(&m_d_nLargeRigidBody, sizeof(int)) );
100:         CudaSafeCall( cudaMalloc(&m_d_largeRigidIndices, m_nLargeRigidBody * sizeof(int)) ); 97:         CudaSafeCall( cudaMalloc(&m_d_largeRigidIndices, m_nLargeRigidBody * sizeof(int)) );
101:  98: 
102:         CudaSafeCall( cudaMalloc(&m_d_ljAddRep, m_nAddTarget * m_nAddTarget * sizeof(double)) ); 
103:         CudaSafeCall( cudaMalloc(&m_d_ljAddAtt, m_nAddTarget * m_nAddTarget * sizeof(double)) ); 
104:  
105:         int counter = 0; 99:         int counter = 0;
106:         while (counter < m_nLargeRigidBody) {100:         while (counter < m_nLargeRigidBody) {
107:                 for (int i = 0; i < m_nRigidBody; ++i) {101:                 for (int i = 0; i < m_nRigidBody; ++i) {
108:                         if (m_nRigidSitesPerBody[i] > 32) {102:                         if (m_nRigidSitesPerBody[i] > 32) {
109:                                 m_largeRigidIndices[counter] = i;103:                                 m_largeRigidIndices[counter] = i;
110:                                 ++counter;104:                                 ++counter;
111:                         }105:                         }
112:                 }106:                 }
113:         }107:         }
114: 108: 
154: 148: 
155:         CudaSafeCall( cudaMemcpy(m_d_isAtomFrozen, m_isAtomFrozen, numDimensions/3 * sizeof(bool), cudaMemcpyHostToDevice) );149:         CudaSafeCall( cudaMemcpy(m_d_isAtomFrozen, m_isAtomFrozen, numDimensions/3 * sizeof(bool), cudaMemcpyHostToDevice) );
156: 150: 
157:         CudaSafeCall( cudaMemcpy(m_d_zerosx, zerosx, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );151:         CudaSafeCall( cudaMemcpy(m_d_zerosx, zerosx, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
158:         CudaSafeCall( cudaMemcpy(m_d_zerosy, zerosy, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );152:         CudaSafeCall( cudaMemcpy(m_d_zerosy, zerosy, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
159:         CudaSafeCall( cudaMemcpy(m_d_zerosz, zerosz, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );153:         CudaSafeCall( cudaMemcpy(m_d_zerosz, zerosz, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
160: 154: 
161:         CudaSafeCall( cudaMemcpy(m_d_nLargeRigidBody, &m_nLargeRigidBody, sizeof(int), cudaMemcpyHostToDevice) );155:         CudaSafeCall( cudaMemcpy(m_d_nLargeRigidBody, &m_nLargeRigidBody, sizeof(int), cudaMemcpyHostToDevice) );
162:         CudaSafeCall( cudaMemcpy(m_d_largeRigidIndices, m_largeRigidIndices, m_nLargeRigidBody * sizeof(int), cudaMemcpyHostToDevice) );156:         CudaSafeCall( cudaMemcpy(m_d_largeRigidIndices, m_largeRigidIndices, m_nLargeRigidBody * sizeof(int), cudaMemcpyHostToDevice) );
163: 157: 
164:         CudaSafeCall( cudaMemcpy(m_d_ljAddRep, m_ljAddRep, m_nAddTarget * m_nAddTarget * sizeof(double), cudaMemcpyHostToDevice) ); 
165:         CudaSafeCall( cudaMemcpy(m_d_ljAddAtt, m_ljAddAtt, m_nAddTarget * m_nAddTarget * sizeof(double), cudaMemcpyHostToDevice) ); 
166:  
167:         delete [] zerosx;158:         delete [] zerosx;
168:         delete [] zerosy;159:         delete [] zerosy;
169:         delete [] zerosz;160:         delete [] zerosz;
170: }161: }
171: 162: 
172: 163: 
173: 164: 
174: CostFunction::~CostFunction() 165: CostFunction::~CostFunction() 
175: {166: {
176:         CudaSafeCall( cudaFree(m_d_zerosx) );167:         CudaSafeCall( cudaFree(m_d_zerosx) );
197:                 CudaSafeCall( cudaFree(m_d_coords) );188:                 CudaSafeCall( cudaFree(m_d_coords) );
198:                 CudaSafeCall( cudaFree(m_d_energy) );189:                 CudaSafeCall( cudaFree(m_d_energy) );
199:                 CudaSafeCall( cudaFree(m_d_gradient) );190:                 CudaSafeCall( cudaFree(m_d_gradient) );
200:         }191:         }
201: 192: 
202:         CudaSafeCall( cudaFree(m_d_isAtomFrozen) );193:         CudaSafeCall( cudaFree(m_d_isAtomFrozen) );
203: 194: 
204:         CudaSafeCall( cudaFree(m_d_nLargeRigidBody) );195:         CudaSafeCall( cudaFree(m_d_nLargeRigidBody) );
205:         CudaSafeCall( cudaFree(m_d_largeRigidIndices) );196:         CudaSafeCall( cudaFree(m_d_largeRigidIndices) );
206: 197: 
207:         CudaSafeCall( cudaFree(m_d_ljAddRep) ); 
208:         CudaSafeCall( cudaFree(m_d_ljAddAtt) ); 
209:  
210:         delete [] m_largeRigidIndices;198:         delete [] m_largeRigidIndices;
211: }199: }
212: 200: 
213: 201: 
214: 202: 
215: // Get/set functions. 203: // Get/set functions. 
216: 204: 
217: size_t CostFunction::getNumDimensions() const205: size_t CostFunction::getNumDimensions() const
218: {206: {
219:         return m_numDimensions;207:         return m_numDimensions;


r32161/cost_function.h 2017-03-22 12:30:28.424245236 +0000 r32160/cost_function.h 2017-03-22 12:30:31.956291835 +0000
 24: // GPU based cost functions can directly inherit from this class 24: // GPU based cost functions can directly inherit from this class
 25: // and implement computeEnergyAndGradient(), the passed pointers are supposed 25: // and implement computeEnergyAndGradient(), the passed pointers are supposed
 26: // to reside in device memory. 26: // to reside in device memory.
 27: class CostFunction 27: class CostFunction
 28: { 28: {
 29:         public: 29:         public:
 30:                 CostFunction(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom,  30:                 CostFunction(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom, 
 31:                                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody,  31:                                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody, 
 32:                                 int *rigidSingles, double *rigidInverse, double *coords, int nSecDiag, bool isAtomisticNotRigid,  32:                                 int *rigidSingles, double *rigidInverse, double *coords, int nSecDiag, bool isAtomisticNotRigid, 
 33:                                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze,  33:                                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze, 
 34:                                 bool isAaConvergence, int nAddTarget, double *ljAddRep, double *ljAddAtt); 34:                                 bool isAaConvergence);
 35:  35: 
 36:                 ~CostFunction(); 36:                 ~CostFunction();
 37:  37: 
 38:                 // Calculates secdiag if R-R minimization and calls rigid body transformations and computeEnergyAndGradient.  38:                 // Calculates secdiag if R-R minimization and calls rigid body transformations and computeEnergyAndGradient. 
 39:                 void basePotential(double *d_x, double *d_f, double *d_gradf); 39:                 void basePotential(double *d_x, double *d_f, double *d_gradf);
 40:  40: 
 41:                 // Implement this method computing both function value d_f 41:                 // Implement this method computing both function value d_f
 42:                 // and gradient d_gradf of your cost function at point d_x. 42:                 // and gradient d_gradf of your cost function at point d_x.
 43:                 virtual void computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf) = 0; 43:                 virtual void computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf) = 0;
 44:  44: 
 66:                 // More accurate calculation of RMS force for rigid body framework (RBF).  66:                 // More accurate calculation of RMS force for rigid body framework (RBF). 
 67:                 void aaConvergence(const double *d_gradf, double *outRms); 67:                 void aaConvergence(const double *d_gradf, double *outRms);
 68:  68: 
 69:                 // Orthogonalise eigenvector to overall rotations and translations.  69:                 // Orthogonalise eigenvector to overall rotations and translations. 
 70:                 void orthogopt(double *d_x, const bool shouldNormalise); 70:                 void orthogopt(double *d_x, const bool shouldNormalise);
 71:  71: 
 72:                 // Freeze specified atoms by setting gradient components to zero.  72:                 // Freeze specified atoms by setting gradient components to zero. 
 73:                 void freeze(double *d_gradf); 73:                 void freeze(double *d_gradf);
 74:  74: 
 75:         protected: 75:         protected:
 76:                 Cublas &m_cublas; 
 77:  
 78:                 const size_t m_numDimensions; // Number of dimensions of function - usually 3*natoms.  76:                 const size_t m_numDimensions; // Number of dimensions of function - usually 3*natoms. 
 79:                 bool m_isColdFusion; // True if cold fusion has occurred.  77:                 bool m_isColdFusion; // True if cold fusion has occurred. 
 80:                 const double m_coldFusionLim; // Energy limit below which we say cold fusion has occurred.  78:                 const double m_coldFusionLim; // Energy limit below which we say cold fusion has occurred. 
 81:                 const int m_nAddTarget; // Target cluster size for addressability.  
 82:                 const double *m_ljAddRep; // Repulsive epsilon matrix for addressability.  
 83:                 const double *m_ljAddAtt; // Attractive epsilon matrix for addressability.  
 84:  
 85:                 int *m_d_nAddTarget; // Target cluster size for addressability, on GPU. 
 86:                 double *m_d_ljAddRep; // Repulsive epsilon matrix for addressability, on GPU. 
 87:                 double *m_d_ljAddAtt; // Attractive epsilon matrix for addressability, on GPU. 
 88:  
 89:  79: 
 90:         private: 80:         private:
 91:                 Printing &m_debugPrinting; 81:                 Printing &m_debugPrinting;
 92:  82: 
 93:                 Timer &m_timer_potential; 83:                 Timer &m_timer_potential;
 94:  84: 
  85:                 Cublas &m_cublas;
  86: 
 95:                 const double  m_aaConvThreshold;      // Threshold under which aaConvergence function is called in RMS force calc. 87:                 const double  m_aaConvThreshold;      // Threshold under which aaConvergence function is called in RMS force calc.
 96:  88: 
 97:                 const double *m_coords;               // The initial coordinates in Bfgsts on CPU, used for R-R minimization. 89:                 const double *m_coords;               // The initial coordinates in Bfgsts on CPU, used for R-R minimization.
 98:  90: 
 99:                 const int     m_nSecDiag;             // Method used to approximate eigenvalue in R-R LBFGS. 91:                 const int     m_nSecDiag;             // Method used to approximate eigenvalue in R-R LBFGS.
100:                 const int     m_nFreeze;              // No. of frozen atoms. 92:                 const int     m_nFreeze;              // No. of frozen atoms.
101:  93: 
102:                 int           m_nCalls;               // No. of potential calls made during this invocation of the CUDA code. 94:                 int           m_nCalls;               // No. of potential calls made during this invocation of the CUDA code.
103:  95: 
104:                 const bool    m_shouldFreeze;         // If true, freeze some specified coordinates. 96:                 const bool    m_shouldFreeze;         // If true, freeze some specified coordinates.
105:  97: 
106:                 const bool   *m_isAtomFrozen;         // Logical array specifying frozen atoms. 98:                 const bool   *m_isAtomFrozen;         // Logical array specifying frozen atoms.
107:  99: 
108:                 bool          m_isBfgsts;             // True if costFunction is being used in a Bfgsts calculation - stops aaConvergence call.100:                 bool          m_isBfgsts;             // True if costFunction is being used in a Bfgsts calculation - stops aaConvergence call.
109:                 bool          m_isRayleighRitz;       // If true, curvature of PES calculated along direction vec at point coords (R-R min).101:                 bool          m_isRayleighRitz;       // If true, curvature of PES calculated along direction vec at point coords (R-R min).
110: 102: 
111:                 const bool    m_isAtomisticNotRigid;  // If false, use rigid body coordinates.103:                 const bool    m_isAtomisticNotRigid;  // If false, use rigid body coordinates.
112:                 const bool    m_isAaConvergence;      // RBF: if true, use more accurate method of calculating RMS force.104:                 const bool    m_isAaConvergence;      // RBF: if true, use more accurate method of calculating RMS force.
113: 105: 
114:                 const double *m_sitesRigidBody;       // RBF: coordinates of the rigid body sites.106:                 const double *m_sitesRigidBody;       // RBF: coordinates of the rigid body sites.
115:                 const double *m_rigidInverse;         // RBF: inverse eigenvalues of the unweighted tensor of gyration. 107:                 const double *m_rigidInverse;         // RBF: inverse eigenvalues of the unweighted tensor of gyration. 
116: 108: 
117:                 const int     m_nDegFreedom;          // RBF: no. of degrees of freedom.109:                 const int     m_nDegFreedom;          // RBF: no. of degrees of freedom.
118:                 const int     m_nRigidBody;           // RBF: no. of rigid bodies.110:                 const int     m_nRigidBody;           // RBF: no. of rigid bodies.
119:                 const int     m_rigidMaxSite;         // RBF: max. no. of sites in a rigid body.111:                 const int     m_rigidMaxSite;         // RBF: max. no. of sites in a rigid body.
120: 112: 
121:                 const int    *m_nRigidSitesPerBody;   // RBF: no. of rigid body sites.113:                 const int    *m_nRigidSitesPerBody;   // RBF: no. of rigid body sites.
122:                 const int    *m_rigidGroups;          // RBF: list of atoms in rigid bodies.114:                 const int    *m_rigidGroups;          // RBF: list of atoms in rigid bodies.
123:                 const int    *m_rigidSingles;         // RBF: list of atoms not in rigid bodies.115:                 const int    *m_rigidSingles;         // RBF: list of atoms not in rigid bodies.
124: 116: 
125:                 int           m_nLargeRigidBody;      // RBF: no. of rigid bodies with more than 32 sites.117:                 int           m_nLargeRigidBody;      // RBF: no. of rigid bodies with more than 32 sites.
126: 118: 
127:                 int          *m_largeRigidIndices;    // RBF: array containing indices for large bodies in nRigidSitesPerBody array. 119:                 int          *m_largeRigidIndices;    // RBF: array containing indices for large bodies in nRigidSitesPerBody array. 
128: 120: 
129: 121: 
130:                 double       *m_d_coords;             // The current coordinates in Bfgsts on GPU, used in R-R minimization. 122:                 double       *m_d_coords;             // The current coordinates in Bfgsts on GPU, used in R-R minimization. 
131:                 double       *m_d_energy;             // The energy at the current coordinates in Bfgsts on GPU, used in R-R minimization. 123:                 double       *m_d_energy;             // The energy at the current coordinates in Bfgsts on GPU, used in R-R minimization. 
132:                 double       *m_d_gradient;           // The gradient at the current coordinates in Bfgsts on GPU, used in R-R minimization. 124:                 double       *m_d_gradient;           // The gradient at the current coordinates in Bfgsts on GPU, used in R-R minimization. 
133:                 double       *m_d_zerosx;             // Array used in cuBLAS operations with ones at all the x positions. 125:                 double       *m_d_zerosx;             // Array used in cuBLAS operations with ones at all the x positions. 
134:                 double       *m_d_zerosy;             // Array used in cuBLAS operations with ones at all the y positions.126:                 double       *m_d_zerosy;             // Array used in cuBLAS operations with ones at all the y positions.
135:                 double       *m_d_zerosz;             // Array used in cuBLAS operations with ones at all the z positions.127:                 double       *m_d_zerosz;             // Array used in cuBLAS operations with ones at all the z positions.
140:                 double       *m_d_rigidInverse;       // RBF: IINVERSE from genrigid, used in aaconvergence subroutine, on GPU.132:                 double       *m_d_rigidInverse;       // RBF: IINVERSE from genrigid, used in aaconvergence subroutine, on GPU.
141:                 double       *m_d_xRigid;             // RBF: rigid body coordinates, on GPU.133:                 double       *m_d_xRigid;             // RBF: rigid body coordinates, on GPU.
142:                 double       *m_d_gkRigid;            // RBF: rigid body gradient, on GPU.134:                 double       *m_d_gkRigid;            // RBF: rigid body gradient, on GPU.
143:                 double       *m_d_gkAtoms;            // RBF: copy of atomistic gradient for use with aaConvergence, on GPU.135:                 double       *m_d_gkAtoms;            // RBF: copy of atomistic gradient for use with aaConvergence, on GPU.
144: 136: 
145:                 int          *m_d_nRigidSitesPerBody; // RBF: no. of rigid body sites, on GPU. 137:                 int          *m_d_nRigidSitesPerBody; // RBF: no. of rigid body sites, on GPU. 
146:                 int          *m_d_rigidGroups;        // RBF: list of atoms in rigid bodies, on GPU. 138:                 int          *m_d_rigidGroups;        // RBF: list of atoms in rigid bodies, on GPU. 
147:                 int          *m_d_rigidSingles;       // RBF: list of atoms not in rigid bodies, on GPU. 139:                 int          *m_d_rigidSingles;       // RBF: list of atoms not in rigid bodies, on GPU. 
148:                 int          *m_d_nRigidBody;         // RBF: no. of rigid bodies, on GPU. 140:                 int          *m_d_nRigidBody;         // RBF: no. of rigid bodies, on GPU. 
149:                 int          *m_d_rigidMaxSite;       // RBF: maximum number of sites in a rigid body, on GPU. 141:                 int          *m_d_rigidMaxSite;       // RBF: maximum number of sites in a rigid body, on GPU. 
150:                 int          *m_d_largeRigidIndices;  // RBF: array containing indices for large bodies in nRigidSitesPerBody array, on GPU. 142:                 int          *m_d_largeRigidIndices;  // RBF: array containing indices for large bodies in nRigidSitesPerBody array, on GPU. 
151:                 int          *m_d_nLargeRigidBody;    // RBF: no. of rigid bodies with more than 32 sites, on GPU. 143:                 int          *m_d_nLargeRigidBody;    // RBF: no. of rigid bodies with more than 32 sites, on GPU. 
 144: 
152: 145: 
153:                 // Rigid body gradient transformations either side of energy/gradient call.146:                 // Rigid body gradient transformations either side of energy/gradient call.
154:                 void rigidTransforms(double *d_x, double *d_f, double *d_gradf);147:                 void rigidTransforms(double *d_x, double *d_f, double *d_gradf);
155: };148: };
156: 149: 
157: #endif /* end of include guard: COST_FUNCTION_H */150: #endif /* end of include guard: COST_FUNCTION_H */


r32161/keywords.f 2017-03-22 12:30:30.640274473 +0000 r32160/keywords.f 2017-03-22 12:30:34.188321280 +0000
733:          MIEF_BOX(1:3) = 1.0D9733:          MIEF_BOX(1:3) = 1.0D9
734:          MIEF_RCUT= 1.0D9734:          MIEF_RCUT= 1.0D9
735:          !735:          !
736:          ! UNDOCUMENTED keywords/parameters736:          ! UNDOCUMENTED keywords/parameters
737:          ! 737:          ! 
738:          TWISTT=.FALSE.738:          TWISTT=.FALSE.
739:          LJADDT=.FALSE.739:          LJADDT=.FALSE.
740:          LJADD2T=.FALSE.740:          LJADD2T=.FALSE.
741:          LJADD3T=.FALSE.741:          LJADD3T=.FALSE.
742:          LJADD4T=.FALSE.742:          LJADD4T=.FALSE.
743:          NADDTARGET=1 
744:          PYADDT=.FALSE.743:          PYADDT=.FALSE.
745:          PYADD2T=.FALSE.744:          PYADD2T=.FALSE.
746:          INVERTPT=.FALSE.745:          INVERTPT=.FALSE.
747:          DNEBEFRAC=0.0D0746:          DNEBEFRAC=0.0D0
748:          MORPHT=.FALSE.747:          MORPHT=.FALSE.
749:          GREATCIRCLET=.FALSE.748:          GREATCIRCLET=.FALSE.
750:          MAXTSENERGY=1.0D100749:          MAXTSENERGY=1.0D100
751:          MAXBARRIER=1.0D100750:          MAXBARRIER=1.0D100
752:          MAXMAXBARRIER=1.0D100751:          MAXMAXBARRIER=1.0D100
753:          ReoptimiseEndpoints=.False.752:          ReoptimiseEndpoints=.False.


r32161/lj_potential.cu 2017-03-22 12:30:29.304256847 +0000 r32160/lj_potential.cu 2017-03-22 12:30:32.844303551 +0000
  2:  *  2:  *
  3:  * File lj_potential.cu: Implementation of class LjPotential, inherited from class CostFunction.   3:  * File lj_potential.cu: Implementation of class LjPotential, inherited from class CostFunction. 
  4:  *  4:  *
  5:  **/  5:  **/
  6:   6: 
  7: #include "cost_function.h"  7: #include "cost_function.h"
  8: #include "potential.h"  8: #include "potential.h"
  9:   9: 
 10: namespace gpu_lj 10: namespace gpu_lj
 11: { 11: {
 12:         // Variables on the GPU.  12:         __constant__ size_t m_d_numDimensions;
 13:  13: 
 14:         __constant__ size_t numDimensions; 14:         // Inner core kernel for LJ.  
 15:         __constant__ int nAddTarget; 15:         // Returns the gradient in x, y, z and energy in w. 
 16:  16:         __device__ double4 ljPair(double4 r1, double4 r2, double4 g)        
 17:  17:         {
 18:         // Kernels. 18:                 double3 r;
 19:  19:                 // Calculate the distance vector. 
 20:         __global__ void fullReduce1(const double *in, double *out, double *d_gradf, const int N, const int roundedMaxSite,  20:                 r.x = r2.x - r1.x;
 21:                         const int blocks); 21:                 r.y = r2.y - r1.y;
 22:         __global__ void kernelComputeEneGrad(const double *d_x, const double *m_d_ljAddRep, const double *m_d_ljAddAtt,  22:                 r.z = r2.z - r1.z;
 23:                         double *d_energyContrib, double *d_gradientContrib); 23: 
 24:         __global__ void deviceReduceKernel(const double *in, double* out, const int N); 24:                 // Distance squared. 
 25: } 25:                 double rsq = r.x * r.x + r.y * r.y + r.z * r.z;
 26:  26: 
 27:  27:                 // Ignore pair if distance is too small, this is usually the self-interaction of the atom. 
 28:  28:                 if (rsq < 1.0e-6)
 29: LjPotential::LjPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, 29:                         return g;
 30:                 int nDegFreedom, int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, 30: 
 31:                 int *rigidGroups, double *sitesRigidBody, int *rigidSingles, double *rigidInverse, 31:                 // Calculate 1/r**2, 1/r**6, 1/r**12. 
 32:                 double *coords, int nSecDiag, bool isAtomisticNotRigid, double aaConvThreshold, 32:                 double ir2  = 1.0 / rsq;
 33:                 double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze, bool isAaConvergence, int nAddTarget, 33:                 double ir6  = ir2*ir2*ir2;
 34:                 double *ljAddRep, double *ljAddAtt) 34:                 double ir12 = ir6*ir6;
 35: : CostFunction(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 35: 
 36:                 nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, coords, 36:                 // Calculate the energy. 
 37:                 nSecDiag, isAtomisticNotRigid, aaConvThreshold, coldFusionLim, shouldFreeze, isAtomFrozen, 37:                 double tmp = 2.0 * (ir12 - ir6);
 38:                 nFreeze, isAaConvergence, nAddTarget, ljAddRep, ljAddAtt) 38:                 g.w += tmp;
 39: { 39:                 // Calculate the gradient. 
 40:         CudaSafeCall( cudaMemcpyToSymbol(gpu_lj::numDimensions, &m_numDimensions, sizeof(size_t)) ); 40:                 tmp = 4.0 * (12.0 * ir12 - 6.0 * ir6) * ir2;
 41:         CudaSafeCall( cudaMemcpyToSymbol(gpu_lj::nAddTarget, &m_nAddTarget, sizeof(int)) ); 41:                 g.x += tmp*r.x;
 42: } 42:                 g.y += tmp*r.y;
 43:  43:                 g.z += tmp*r.z;
 44:  
 45:  
 46: void LjPotential::computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf) 
 47: { 
 48:         using namespace gpu_lj; 
 49:  
 50:         dim3 blockDim; 
 51:         dim3 gridDim; 
 52:  
 53:         // Set launch parameters for main potential kernel.  
 54:         int numThreads = (m_numDimensions/3)*(m_numDimensions/3); 
 55:         blockDim.x = 512; 
 56:         int numBlocks = (numThreads + blockDim.x - 1)/blockDim.x; 
 57:         gridDim.x = numBlocks; 
 58:  
 59:         double *d_energyContrib; 
 60:         double *d_gradientContrib; 
 61:  
 62:         CudaSafeCall( cudaMalloc(&d_energyContrib, numThreads * sizeof(double)) ); 
 63:         CudaSafeCall( cudaMalloc(&d_gradientContrib, 3 * numThreads * sizeof(double)) ); 
 64:  
 65:         gpu_lj::kernelComputeEneGrad<<<gridDim, blockDim>>>(d_x, m_d_ljAddRep, m_d_ljAddAtt, d_energyContrib, d_gradientContrib); 
 66:         CudaCheckError(); 
 67:         cudaDeviceSynchronize(); 
 68:  
 69:  
 70:         // Summation of energy contributions 
 71:         int threads = 512; 
 72:         int blocks = min((numThreads + threads - 1) / threads, 1024); 
 73:  
 74:         double *d_out; 
 75:  
 76:         CudaSafeCall( cudaMalloc(&d_out, blocks * sizeof(double)) ); 
 77:  
 78:         gpu_lj::deviceReduceKernel<<<blocks, threads>>>(d_energyContrib, d_out, numThreads); 
 79:         CudaCheckError(); 
 80:         cudaDeviceSynchronize(); 
 81:  
 82:         gpu_lj::deviceReduceKernel<<<1, 1024>>>(d_out, d_out, blocks); 
 83:         CudaCheckError(); 
 84:         cudaDeviceSynchronize(); 
 85:  
 86:         CudaSafeCall( cudaMemcpy(d_f, d_out, sizeof(double), cudaMemcpyDeviceToDevice) ); 
 87:  
 88:  
 89:         // Summation of gradient contributions 
 90:         int blockSize = 1024; 
 91:         blockDim.x = blockSize; 
 92:  
 93:         // Round m_numDimensions to nearest block size so that per block reduction doesn't sum components from different atoms.  
 94:         int roundedMaxSite = ((m_numDimensions/3 + blockSize - 1)/blockSize)*blockSize; 
 95:  
 96:         numThreads = roundedMaxSite*m_numDimensions/3; 
 97:  
 98:         blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
 99:         gridDim.x = blocks; 
100:  
101:         double *d_outArray; 
102:  
103:         CudaSafeCall( cudaMalloc(&d_outArray, 3 * blocks * sizeof(double)) ); 
104:  
105:         gpu_lj::fullReduce1<<<gridDim, blockDim>>>(d_gradientContrib, d_outArray, d_gradf, m_numDimensions/3, roundedMaxSite,  
106:                         blocks); 
107:         CudaCheckError(); 
108:         cudaDeviceSynchronize(); 
109:  
110:         while (blocks > m_numDimensions/3) { 
111:                 int oldBlocks = blocks; 
112:  
113:                 roundedMaxSite = ((blocks/(m_numDimensions/3) + blockSize - 1)/blockSize)*blockSize; 
114:                 numThreads = roundedMaxSite*m_numDimensions/3; 
115:  
116:                 blocks = (numThreads + blockDim.x - 1)/blockDim.x; 
117:                 gridDim.x = blocks; 
118:  44: 
119:                 gpu_lj::fullReduce1<<<gridDim, blockDim>>>(d_outArray, d_outArray, d_gradf, oldBlocks/(m_numDimensions/3),  45:                 return g;
120:                                 roundedMaxSite, blocks); 
121:                 CudaCheckError(); 
122:                 cudaDeviceSynchronize(); 
123:         } 46:         }
124:  47: 
125:         CudaSafeCall( cudaFree(d_energyContrib) ); 48:         // Perform the calculation for one block. 
126:         CudaSafeCall( cudaFree(d_gradientContrib) ); 49:         __device__ double4 calculateTile(double4 r, double4 g)
127:         CudaSafeCall( cudaFree(d_out) ); 50:         {
128:         CudaSafeCall( cudaFree(d_outArray) ); 51:                 extern __shared__ double4 shared_r[];
129: } 52:                 for (size_t i = 0; i < (m_d_numDimensions/3); i++) {
130:  53:                         g=gpu_lj::ljPair(r, shared_r[i], g);
131:  
132:  
133: namespace gpu_lj  
134: { 
135:         // Undocumented CUDA 6.5 function copied and pasted here as not present in CUDA 6.0 and below.  
136:         static __device__ __inline__ 
137:                 double __myshfl_down(double var, unsigned int delta, int width=warpSize) { 
138:                         float lo, hi; 
139:                         asm volatile("mov.b64 {%0,%1}, %2;" : "=f"(lo), "=f"(hi) : "d"(var)); 
140:                         hi = __shfl_down(hi, delta, width); 
141:                         lo = __shfl_down(lo, delta, width); 
142:                         asm volatile("mov.b64 %0, {%1,%2};" : "=d"(var) : "f"(lo), "f"(hi)); 
143:                         return var; 
144:                 } 
145:  
146:         __inline__ __device__ 
147:                 double warpReduceSum(double val) { 
148:                         for (int offset = warpSize/2; offset > 0; offset /= 2)  
149:                                 val += __myshfl_down(val, offset); 
150:                         return val; 
151:                 } 54:                 }
  55:                 return g;
  56:         }
152:  57: 
153:         __inline__ __device__ 58:         __global__ void kernelComputeEneGrad(const double *d_x, double *d_y, double *d_grad)
154:                 double blockReduceSum(double val) { 59:         {
155:                         static __shared__ double shared[32]; // Shared mem for 32 partial sums 60:                 extern __shared__ double4 shared_r[];
156:                         int lane = threadIdx.x % warpSize; 61:                 double4 myPosition; // Each thread deals with one atom
157:                         int wid = threadIdx.x / warpSize; 62:                 double4 g = {0.0, 0.0, 0.0, 0.0};        
158:  
159:                         val = warpReduceSum(val);     // Each warp performs partial reduction 
160:  
161:                         if (lane==0) shared[wid]=val; // Write reduced value to shared memory 
162:  
163:                         __syncthreads();              // Wait for all partial reductions 
164:  
165:                         //read from shared memory only if that warp existed 
166:                         val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0; 
167:  63: 
168:                         if (wid==0) val = warpReduceSum(val); //Final reduce within first warp 64:                 // Current test version of LJ uses only one block for efficient shared memory use. 
  65:                 int tid = threadIdx.x;
169:  66: 
170:                         return val; 67:                 // Make sure w is initialised to 0 on ALL shared memory so reduction is safe. 
171:                 } 68:                 shared_r[tid].w = 0.0;
172:  69: 
173:         __global__ void deviceReduceKernel(const double *in, double* out, const int N) { 70:                 // Thread numbers greater than the number of atoms do nothing. 
174:                 double sum = 0; 71:                 if (tid < (m_d_numDimensions/3)) {
175:                 //reduce multiple elements per thread. 72:                         // Read the coordinates from global memory into memory local to each thread. 
176:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; 73:                         myPosition.x = d_x[3*tid+0];
177:                                 i < N; 74:                         myPosition.y = d_x[3*tid+1];
178:                                 i += blockDim.x * gridDim.x) { 75:                         myPosition.z = d_x[3*tid+2];
179:                         sum += in[i]; 76: 
180:                 } 77:                         myPosition.w = 0.0;
181:                 sum = blockReduceSum(sum); 78: 
182:                 if (threadIdx.x==0) 79:                         // Also copy coordinates from global memory into shared memory. 
183:                         out[blockIdx.x]=sum; 80:                         shared_r[tid].x = myPosition.x; 
184:         } 81:                         shared_r[tid].y = myPosition.y; 
185:  82:                         shared_r[tid].z = myPosition.z; 
186:         __global__ void fullReduce1(const double *in, double* out, double *d_gradf, const int N, const int roundedMaxSite,  
187:                         const int blocks) 
188:         { 
189:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
190:  83: 
191:                 while (tid < (roundedMaxSite*numDimensions/3)) { 84:                         // Wait until copy to shared memory completed, 
  85:                         __syncthreads();
192:  86: 
193:                         double3 elements; 87:                         // Carry out energy and gradient calculations. 
194:                         elements.x = 0.0; 88:                         g = gpu_lj::calculateTile(myPosition, g);
195:                         elements.y = 0.0; 
196:                         elements.z = 0.0; 
197:  
198:                         int index = ((tid + roundedMaxSite) % roundedMaxSite); 
199:                         int thisAtom = tid/roundedMaxSite; 
200:  
201:                         if (index < N) { 
202:                                 int var = index + thisAtom*N; 
203:                                 elements.x = in[3*var+0]; 
204:                                 elements.y = in[3*var+1]; 
205:                                 elements.z = in[3*var+2]; 
206:                         } 
207:  89: 
  90:                         // Wait until the calculations are done. 
208:                         __syncthreads(); 91:                         __syncthreads();
209:  92: 
210:                         elements.x = blockReduceSum(elements.x); 93:                         // Copy energy in w to shared memory for reduction and copy gradient to global memory. 
211:                         elements.y = blockReduceSum(elements.y); 94:                         shared_r[tid].w = g.w;
212:                         elements.z = blockReduceSum(elements.z); 95:                         d_grad[3*tid+0] = g.x;
  96:                         d_grad[3*tid+1] = g.y;
  97:                         d_grad[3*tid+2] = g.z;
213:  98: 
214:                         __syncthreads(); 99:                         __syncthreads();
215: 100: 
216:                         if (threadIdx.x==0) {101:                 }
217:                                 if (blocks == numDimensions/3) {102:                 // Shared memory reduction to add up all energy contributions to get a total energy. 
218:                                         d_gradf[3*thisAtom+0] = elements.x;103:                 int i = blockDim.x/2;
219:                                         d_gradf[3*thisAtom+1] = elements.y;104:                 while (i != 0) {
220:                                         d_gradf[3*thisAtom+2] = elements.z;105:                         if (tid < i) {
221:                                 }106:                                 shared_r[tid].w += shared_r[tid+i].w;
222:                                 else { 
223:                                         out[3*blockIdx.x+0] = elements.x; 
224:                                         out[3*blockIdx.x+1] = elements.y; 
225:                                         out[3*blockIdx.x+2] = elements.z; 
226:                                 } 
227:                         }107:                         }
 108:                         __syncthreads();
 109:                         i /= 2;
 110:                 }
228: 111: 
229:                         tid += blockDim.x * gridDim.x;112:                 if (tid==0) {
 113:                         // Write total energy back to global memory. 
 114:                         *d_y = shared_r[0].w;
230:                 }115:                 }
231:         }116:         }
 117: }
232: 118: 
233:         __global__ void kernelComputeEneGrad(const double *d_x, const double *m_d_ljAddRep, const double *m_d_ljAddAtt,  
234:                         double *d_energyContrib, double *d_gradientContrib) 
235:         { 
236:                 // Each thread deals with one atom interacting with one other atom. 
237:  
238:                 int tid = blockIdx.x * blockDim.x + threadIdx.x; 
239: 119: 
240:                 while (tid < (numDimensions/3)*(numDimensions/3)) { 
241: 120: 
242:                         int refAtom = tid / (numDimensions/3); // Integer division rounds down.121: LjPotential::LjPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, 
243:                         int myAtom = tid % (numDimensions/3);122:                 int nDegFreedom, int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, 
 123:                 int *rigidGroups, double *sitesRigidBody, int *rigidSingles, double *rigidInverse, 
 124:                 double *coords, int nSecDiag, bool isAtomisticNotRigid, double aaConvThreshold, 
 125:                 double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze, bool isAaConvergence)
 126: : CostFunction(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 
 127:                 nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, coords, 
 128:                 nSecDiag, isAtomisticNotRigid, aaConvThreshold, coldFusionLim, shouldFreeze, isAtomFrozen, 
 129:                 nFreeze, isAaConvergence)
 130: {
 131:         // Copy 3*natoms to GPU
 132:         CudaSafeCall( cudaMemcpyToSymbol(gpu_lj::m_d_numDimensions, &m_numDimensions, sizeof(size_t)) );
 133: }
244: 134: 
245:                         // Read the coordinates from global memory into memory local to each thread.  
246:                         double myPositionX = d_x[3*myAtom+0]; 
247:                         double myPositionY = d_x[3*myAtom+1]; 
248:                         double myPositionZ = d_x[3*myAtom+2]; 
249:  
250:                         double refPositionX = d_x[3*refAtom+0]; 
251:                         double refPositionY = d_x[3*refAtom+1]; 
252:                         double refPositionZ = d_x[3*refAtom+2]; 
253: 135: 
254:                         int mj1 = refAtom % nAddTarget; 
255:                         int mj2 = myAtom % nAddTarget; 
256: 136: 
257:                         // Carry out energy and gradient calculations. 137: void LjPotential::computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf)
 138: {
 139:         dim3 blockDim;
 140:         dim3 gridDim;
258: 141: 
259:                         // Calculate the distance vector. 142:         // Set launch parameters for main potential kernel. 
260:                         double3 r;143:         // No. of threads (1024 is max. per block). 
261:                         r.x = myPositionX - refPositionX;144:         blockDim.x = 1024;
262:                         r.y = myPositionY - refPositionY;145:         int numBlocks = ((m_numDimensions/3) + blockDim.x - 1)/blockDim.x;
263:                         r.z = myPositionZ - refPositionZ;146:         // Number of blocks (one). 
264: 147:         gridDim.x = numBlocks;
265:                         // Distance squared.  
266:                         double rsq = r.x * r.x + r.y * r.y + r.z * r.z; 
267:  
268:                         double attractive = m_d_ljAddAtt[mj2+mj1*nAddTarget]; 
269:                         double repulsive = m_d_ljAddRep[mj2+mj1*nAddTarget]; 
270:  
271:                         d_energyContrib[tid] = 0.0; 
272:                         d_gradientContrib[3*tid+0] = 0.0; 
273:                         d_gradientContrib[3*tid+1] = 0.0; 
274:                         d_gradientContrib[3*tid+2] = 0.0; 
275:  
276:                         // Ignore pair if distance is too small, this is usually the self-interaction of the atom.  
277:                         if (rsq >= 1.0e-6) { 
278:                                 // Calculate 1/r**2, 1/r**6, 1/r**12.  
279:                                 double ir2  = 1.0 / rsq; 
280:                                 double ir6  = ir2*ir2*ir2; 
281:                                 double ir12 = ir6*ir6; 
282:  
283:                                 // Calculate the energy.  
284:                                 d_energyContrib[tid] = 2.0 * (ir12*repulsive-ir6*attractive); 
285:                                 // Calculate the gradient.  
286:                                 double gradcontrib = 4.0 * (12.0 * ir12 * repulsive - 
287:                                                 6.0 * ir6 * attractive) * ir2; 
288:                                 d_gradientContrib[3*tid+0] = gradcontrib * r.x; 
289:                                 d_gradientContrib[3*tid+1] = gradcontrib * r.y; 
290:                                 d_gradientContrib[3*tid+2] = gradcontrib * r.z; 
291:                         } 
292: 148: 
293:                         tid += blockDim.x * gridDim.x;149:         int sharedSize;
294:                 }150:         // Size of shared memory. 
 151:         sharedSize = sizeof(double4) * (blockDim.x);
 152: 
 153:         // If shared size bigger than device shared memory, warn and exit. 
 154:         int device;
 155:         cudaGetDevice(&device);
 156: 
 157:         cudaDeviceProp props;
 158:         cudaGetDeviceProperties(&props, device);
 159: 
 160:         if (sharedSize > props.sharedMemPerBlock) {
 161:                 std::cerr << "The amount of shared memory requested exceeds availability. " << std::endl;
 162:                 exit(EXIT_FAILURE);
295:         }163:         }
296: 164: 
 165:         gpu_lj::kernelComputeEneGrad<<<gridDim, blockDim, sharedSize>>>(d_x, d_f, d_gradf);
 166:         CudaCheckError();
 167:         cudaDeviceSynchronize();
297: }168: }


r32161/modcudabfgsts.f90 2017-03-22 12:30:30.856277323 +0000 r32160/modcudabfgsts.f90 2017-03-22 12:30:34.408324183 +0000
  1: MODULE MODCUDABFGSTS  1: MODULE MODCUDABFGSTS
  2:   2: 
  3: USE COMMONS, ONLY : NATOMS, DEBUG, NOPT, STPMAX, REDOPATH, IVEC  3: USE COMMONS, ONLY : NATOMS, DEBUG, NOPT, STPMAX, REDOPATH, IVEC
  4:   4: 
  5: USE KEY, ONLY : CUDAPOT, GMAX, CONVU, CONVR, MINMAX, MXSTP, MAXXBFGS, XMAXERISE, CUDATIMET, CFUSIONT, COLDFUSIONLIMIT, XDGUESS, &   5: USE KEY, ONLY : CUDAPOT, GMAX, CONVU, CONVR, MINMAX, MXSTP, MAXXBFGS, XMAXERISE, CUDATIMET, CFUSIONT, COLDFUSIONLIMIT, XDGUESS, & 
  6:                 XMUPDATE, HINDEX, MASST, SHIFTED, CHECKINDEX, REVERSEUPHILLT, GRADSQ, FREEZE, READV, DUMPV, FIXD, EIGENONLY, &   6:                 XMUPDATE, HINDEX, MASST, SHIFTED, CHECKINDEX, REVERSEUPHILLT, GRADSQ, FREEZE, READV, DUMPV, FIXD, EIGENONLY, & 
  7:                 NSECDIAG, EVPC, CEIG, NEVS, CHECKOVERLAPT, CHECKNEGATIVET, REOPT, BFGSSTEP, PUSHCUT, PUSHOFF, MAXMAX, TRAD, &   7:                 NSECDIAG, EVPC, CEIG, NEVS, CHECKOVERLAPT, CHECKNEGATIVET, REOPT, BFGSSTEP, PUSHCUT, PUSHOFF, MAXMAX, TRAD, & 
  8:                 NBFGSMAX1, NBFGSMAX2, BFGSTSTOL, REPELTST, MUPDATE, MAXBFGS, MAXERISE, DGUESS, NOIT, VARIABLES, NOHESS, FROZEN, &   8:                 NBFGSMAX1, NBFGSMAX2, BFGSTSTOL, REPELTST, MUPDATE, MAXBFGS, MAXERISE, DGUESS, NOIT, VARIABLES, NOHESS, FROZEN, & 
  9:                 NFREEZE, PV, FIXAFTER, CPMD, PRINTPTS, TWOENDS, DUMPMAG, NSTEPMIN, AMBER12T, LJADD3T, NADDTARGET, LJADDREP, LJADDATT  9:                 NFREEZE, PV, FIXAFTER, CPMD, PRINTPTS, TWOENDS, DUMPMAG, NSTEPMIN, AMBER12T
 10:  10: 
 11: USE GENRIGID, ONLY : RIGIDINIT, ATOMRIGIDCOORDT, DEGFREEDOMS, NRIGIDBODY, NSITEPERBODY, RIGIDGROUPS, MAXSITE, SITESRIGIDBODY, & 11: USE GENRIGID, ONLY : RIGIDINIT, ATOMRIGIDCOORDT, DEGFREEDOMS, NRIGIDBODY, NSITEPERBODY, RIGIDGROUPS, MAXSITE, SITESRIGIDBODY, &
 12:                      RIGIDSINGLES, IINVERSE, AACONVERGENCET 12:                      RIGIDSINGLES, IINVERSE, AACONVERGENCET
 13:  13: 
 14: USE ZWK, ONLY : NUP, ZWORK 14: USE ZWK, ONLY : NUP, ZWORK
 15:  15: 
 16: USE MODNEB, ONLY : NEWCONNECTT, NEWNEBT 16: USE MODNEB, ONLY : NEWCONNECTT, NEWNEBT
 17:  17: 
 18: USE INTCOMMONS, ONLY : BONDSFROMFILE 18: USE INTCOMMONS, ONLY : BONDSFROMFILE
 19:  19: 
 25:  25: 
 26: IMPLICIT NONE 26: IMPLICIT NONE
 27:  27: 
 28: INTERFACE 28: INTERFACE
 29:     SUBROUTINE CUDA_BFGSTS(C_NATOMS, C_COORDS, C_CEIG, C_MFLAG, C_ENERGY, C_NEVS, ITMAX, C_ITER, C_MAXXBFGS, C_XMAXERISE, C_RMS, & 29:     SUBROUTINE CUDA_BFGSTS(C_NATOMS, C_COORDS, C_CEIG, C_MFLAG, C_ENERGY, C_NEVS, ITMAX, C_ITER, C_MAXXBFGS, C_XMAXERISE, C_RMS, &
 30:                           C_CUDAPOT, C_DEBUG, C_CUDATIMET, NCALLS, C_CFUSIONT, C_COLDFUSIONLIMIT, C_XDGUESS, C_XMUPDATE, & 30:                           C_CUDAPOT, C_DEBUG, C_CUDATIMET, NCALLS, C_CFUSIONT, C_COLDFUSIONLIMIT, C_XDGUESS, C_XMUPDATE, &
 31:                           C_EVALMIN, C_VECS, C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, & 31:                           C_EVALMIN, C_VECS, C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, &
 32:                           C_MAXSITE, C_SITESRIGIDBODY, C_RIGIDSINGLES, C_IINVERSE, C_NSECDIAG, C_EVPC, & 32:                           C_MAXSITE, C_SITESRIGIDBODY, C_RIGIDSINGLES, C_IINVERSE, C_NSECDIAG, C_EVPC, &
 33:                           C_PUSHCUT, C_PUSHOFF, C_STPMAX, C_MAXMAX, C_MINMAX, C_TRAD, C_NBFGSMAX1, C_NBFGSMAX2, C_BFGSTSTOL, & 33:                           C_PUSHCUT, C_PUSHOFF, C_STPMAX, C_MAXMAX, C_MINMAX, C_TRAD, C_NBFGSMAX1, C_NBFGSMAX2, C_BFGSTSTOL, &
 34:                           C_MUPDATE, C_GMAX, C_MAXBFGS, C_MAXERISE, C_DGUESS, C_CONVU, C_FREEZE, C_FROZEN, &  34:                           C_MUPDATE, C_GMAX, C_MAXBFGS, C_MAXERISE, C_DGUESS, C_CONVU, C_FREEZE, C_FROZEN, & 
 35:                           C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET, C_NADDTARGET, C_LJADDREP, &  35:                           C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET)BIND(C,NAME="setup_bfgsts")
 36:                           C_LJADDATT)BIND(C,NAME="setup_bfgsts") 
 37:  36: 
 38:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR 37:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR
 39:  38: 
 40:         INTEGER(KIND=C_INT), INTENT(IN) :: C_NATOMS, & ! No. of atoms 39:         INTEGER(KIND=C_INT), INTENT(IN) :: C_NATOMS, & ! No. of atoms
 41:                                            ITMAX, & ! Max. no. of BFGSTS iterations allowed 40:                                            ITMAX, & ! Max. no. of BFGSTS iterations allowed
 42:                                            C_DEGFREEDOMS, & ! Rigid Body Framework (RBF): no. of degrees of freedom 41:                                            C_DEGFREEDOMS, & ! Rigid Body Framework (RBF): no. of degrees of freedom
 43:                                            C_NRIGIDBODY, & ! RBF: no. of rigid bodies 42:                                            C_NRIGIDBODY, & ! RBF: no. of rigid bodies
 44:                                            C_MAXSITE, & ! RBF: max. no. of sites in a rigid body 43:                                            C_MAXSITE, & ! RBF: max. no. of sites in a rigid body
 45:                                            C_XMUPDATE, & ! History size for Rayleigh-Ritz (R-R) LBFGS minimization 44:                                            C_XMUPDATE, & ! History size for Rayleigh-Ritz (R-R) LBFGS minimization
 46:                                            C_NSECDIAG, & ! Method used to approximate eigenvalue in R-R LBFGS 45:                                            C_NSECDIAG, & ! Method used to approximate eigenvalue in R-R LBFGS
 47:                                            C_NEVS, & ! Max. no. of steps allowed in R-R LBFGS  46:                                            C_NEVS, & ! Max. no. of steps allowed in R-R LBFGS 
 48:                                            C_NBFGSMAX1, & ! No. of LBFGS steps allowed in restricted minimization (grad projected) 47:                                            C_NBFGSMAX1, & ! No. of LBFGS steps allowed in restricted minimization (grad projected)
 49:                                            C_NBFGSMAX2, & ! No. of LBFGS steps allowed in tangent space after eigenvalue converged 48:                                            C_NBFGSMAX2, & ! No. of LBFGS steps allowed in tangent space after eigenvalue converged
 50:                                            C_MUPDATE, & ! History size for restricted LBFGS minimization 49:                                            C_MUPDATE, & ! History size for restricted LBFGS minimization
 51:                                            C_NFREEZE, & ! No. of frozen atoms 50:                                            C_NFREEZE ! No. of frozen atoms
 52:                                            C_NADDTARGET ! Target cluster size (addressability) 
 53:  51: 
 54:         INTEGER(KIND=C_INT), DIMENSION(C_NRIGIDBODY), INTENT(IN) :: C_NSITEPERBODY ! RBF: no. of rigid body sites 52:         INTEGER(KIND=C_INT), DIMENSION(C_NRIGIDBODY), INTENT(IN) :: C_NSITEPERBODY ! RBF: no. of rigid body sites
 55:  53: 
 56:         INTEGER(KIND=C_INT), DIMENSION(C_MAXSITE*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDGROUPS ! RBF: list of atoms in rigid bodies 54:         INTEGER(KIND=C_INT), DIMENSION(C_MAXSITE*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDGROUPS ! RBF: list of atoms in rigid bodies
 57:  55: 
 58:         INTEGER(KIND=C_INT), DIMENSION(C_DEGFREEDOMS/3 - 2*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDSINGLES ! RBF: list of atoms not in rigid bodies 56:         INTEGER(KIND=C_INT), DIMENSION(C_DEGFREEDOMS/3 - 2*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDSINGLES ! RBF: list of atoms not in rigid bodies
 59:  57: 
 60:         INTEGER(KIND=C_INT), INTENT(OUT) :: C_ITER, & ! No. of BFGSTS iterations done 58:         INTEGER(KIND=C_INT), INTENT(OUT) :: C_ITER, & ! No. of BFGSTS iterations done
 61:                                             NCALLS ! No. of potential calls made during this call to BFGSTS  59:                                             NCALLS ! No. of potential calls made during this call to BFGSTS 
 62:  60: 
 76:                                            C_GMAX, & ! Convergence tolerance for RMS force for restricted LBFGS minimization 74:                                            C_GMAX, & ! Convergence tolerance for RMS force for restricted LBFGS minimization
 77:                                            C_MAXBFGS, & ! Max. step size allowed for restricted LBFGS 75:                                            C_MAXBFGS, & ! Max. step size allowed for restricted LBFGS
 78:                                            C_MAXERISE, & ! Max. energy rise allowed for restricted LBFGS 76:                                            C_MAXERISE, & ! Max. energy rise allowed for restricted LBFGS
 79:                                            C_DGUESS, & ! Initial guess for restricted LBFGS inverse Hessian diagonal elements 77:                                            C_DGUESS, & ! Initial guess for restricted LBFGS inverse Hessian diagonal elements
 80:                                            C_CONVU ! Can only converge if eigenvector following step length magnitude falls below this value 78:                                            C_CONVU ! Can only converge if eigenvector following step length magnitude falls below this value
 81:  79: 
 82:         REAL(KIND=C_DOUBLE), DIMENSION(C_MAXSITE*3*C_NRIGIDBODY), INTENT(IN) :: C_SITESRIGIDBODY ! RBF: coordinates of the rigid body sites 80:         REAL(KIND=C_DOUBLE), DIMENSION(C_MAXSITE*3*C_NRIGIDBODY), INTENT(IN) :: C_SITESRIGIDBODY ! RBF: coordinates of the rigid body sites
 83:  81: 
 84:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration 82:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration
 85:  83: 
 86:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 
 87:  
 88:         REAL(KIND=C_DOUBLE), DIMENSION(3*C_NATOMS), INTENT(INOUT) :: C_COORDS, & ! Coordinates 84:         REAL(KIND=C_DOUBLE), DIMENSION(3*C_NATOMS), INTENT(INOUT) :: C_COORDS, & ! Coordinates
 89:                                                                      C_VECS ! Lowest eigenvector 85:                                                                      C_VECS ! Lowest eigenvector
 90:  86: 
 91:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy at coordinates C_COORDS 87:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy at coordinates C_COORDS
 92:                                             C_RMS, & ! RMS force 88:                                             C_RMS, & ! RMS force
 93:                                             C_EVALMIN, & ! Smallest eigenvalue 89:                                             C_EVALMIN, & ! Smallest eigenvalue
 94:                                             POTENTIALTIME ! Time taken in calculating potential 90:                                             POTENTIALTIME ! Time taken in calculating potential
 95:  91: 
 96:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info. 92:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info.
 97:                                             C_CUDATIMET, & ! If true, print timing info. 93:                                             C_CUDATIMET, & ! If true, print timing info.
108: 104: 
109:     END SUBROUTINE CUDA_BFGSTS105:     END SUBROUTINE CUDA_BFGSTS
110: END INTERFACE106: END INTERFACE
111:  107:  
112: CONTAINS108: CONTAINS
113: 109: 
114:     SUBROUTINE CUDA_BFGSTS_WRAPPER(ITMAX,COORDS,ENERGY,MFLAG,RMS,EVALMIN,VECS,ITER)110:     SUBROUTINE CUDA_BFGSTS_WRAPPER(ITMAX,COORDS,ENERGY,MFLAG,RMS,EVALMIN,VECS,ITER)
115: 111: 
116:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_BFGSTS are converted directly112:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_BFGSTS are converted directly
117:         INTEGER(KIND=C_INT) :: C_NATOMS, ITMAX, C_ITER, NCALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, C_XMUPDATE, &113:         INTEGER(KIND=C_INT) :: C_NATOMS, ITMAX, C_ITER, NCALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, C_XMUPDATE, &
118:                                C_NSECDIAG, C_NEVS, C_NBFGSMAX1, C_NBFGSMAX2, C_MUPDATE, C_NFREEZE, C_NADDTARGET114:                                C_NSECDIAG, C_NEVS, C_NBFGSMAX1, C_NBFGSMAX2, C_MUPDATE, C_NFREEZE
119:         INTEGER(KIND=C_INT), DIMENSION(NRIGIDBODY) :: C_NSITEPERBODY115:         INTEGER(KIND=C_INT), DIMENSION(NRIGIDBODY) :: C_NSITEPERBODY
120:         INTEGER(KIND=C_INT), DIMENSION(MAXSITE*NRIGIDBODY) :: C_RIGIDGROUPS116:         INTEGER(KIND=C_INT), DIMENSION(MAXSITE*NRIGIDBODY) :: C_RIGIDGROUPS
121:         INTEGER(KIND=C_INT), DIMENSION(DEGFREEDOMS/3 - 2*NRIGIDBODY) :: C_RIGIDSINGLES117:         INTEGER(KIND=C_INT), DIMENSION(DEGFREEDOMS/3 - 2*NRIGIDBODY) :: C_RIGIDSINGLES
122: 118: 
123:         REAL(KIND=C_DOUBLE) :: C_CEIG, C_ENERGY, C_RMS, C_XDGUESS, C_EVALMIN, C_MAXXBFGS, C_XMAXERISE, C_COLDFUSIONLIMIT, &119:         REAL(KIND=C_DOUBLE) :: C_CEIG, C_ENERGY, C_RMS, C_XDGUESS, C_EVALMIN, C_MAXXBFGS, C_XMAXERISE, C_COLDFUSIONLIMIT, &
124:                                C_EVPC, C_PUSHCUT, C_PUSHOFF, C_STPMAX, C_MAXMAX, C_MINMAX, C_TRAD, C_BFGSTSTOL, C_GMAX, &120:                                C_EVPC, C_PUSHCUT, C_PUSHOFF, C_STPMAX, C_MAXMAX, C_MINMAX, C_TRAD, C_BFGSTSTOL, C_GMAX, &
125:                                C_MAXBFGS, C_MAXERISE, C_DGUESS, C_CONVU, POTENTIALTIME121:                                C_MAXBFGS, C_MAXERISE, C_DGUESS, C_CONVU, POTENTIALTIME
126:         REAL(KIND=C_DOUBLE), DIMENSION(MAXSITE*3*NRIGIDBODY) :: C_SITESRIGIDBODY 122:         REAL(KIND=C_DOUBLE), DIMENSION(MAXSITE*3*NRIGIDBODY) :: C_SITESRIGIDBODY 
127:         REAL(KIND=C_DOUBLE), DIMENSION(NRIGIDBODY*3*3) :: C_IINVERSE123:         REAL(KIND=C_DOUBLE), DIMENSION(NRIGIDBODY*3*3) :: C_IINVERSE
128:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT 
129:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: C_COORDS, C_VECS124:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: C_COORDS, C_VECS
130: 125: 
131:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_CFUSIONT, C_ATOMRIGIDCOORDT, C_FREEZE, C_AACONVERGENCET126:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_CFUSIONT, C_ATOMRIGIDCOORDT, C_FREEZE, C_AACONVERGENCET
132:         LOGICAL(KIND=C_BOOL), DIMENSION(NATOMS) :: C_FROZEN127:         LOGICAL(KIND=C_BOOL), DIMENSION(NATOMS) :: C_FROZEN
133: 128: 
134:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT129:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
135: 130: 
136:         ! Variables passed as *arguments through this wrapper* (not common) with intent out for CUDA_BFGSTS are not passed into it131:         ! Variables passed as *arguments through this wrapper* (not common) with intent out for CUDA_BFGSTS are not passed into it
137:         ! Therefore uninitialised C types are passed in and converted types are copied back after the call132:         ! Therefore uninitialised C types are passed in and converted types are copied back after the call
138: 133: 
239:         END DO234:         END DO
240: 235: 
241:         DO K = 1,3236:         DO K = 1,3
242:             DO J = 1,3237:             DO J = 1,3
243:                 DO I = 1,NRIGIDBODY238:                 DO I = 1,NRIGIDBODY
244:                    C_IINVERSE((K - 1)*3*NRIGIDBODY + (J - 1)*NRIGIDBODY + I) = IINVERSE(I,J,K)239:                    C_IINVERSE((K - 1)*3*NRIGIDBODY + (J - 1)*NRIGIDBODY + I) = IINVERSE(I,J,K)
245:                 END DO240:                 END DO
246:             END DO241:             END DO
247:         END DO242:         END DO
248: 243: 
249:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN 
250:             DO J = 1,NADDTARGET 
251:                 DO I = 1,NADDTARGET 
252:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J) 
253:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J) 
254:                 END DO 
255:             END DO 
256:         ELSE 
257:             C_LJADDREP(:) = 1.0 
258:             C_LJADDATT(:) = 1.0 
259:         END IF 
260:  
261:         DO I = 1,3*NATOMS244:         DO I = 1,3*NATOMS
262:             C_COORDS(I) = COORDS(I)245:             C_COORDS(I) = COORDS(I)
263:         END DO246:         END DO
264: 247: 
265:         DO I = 1,NATOMS248:         DO I = 1,NATOMS
266:             C_FROZEN(I) = FROZEN(I)249:             C_FROZEN(I) = FROZEN(I)
267:         END DO250:         END DO
268: 251: 
269:         C_NATOMS = NATOMS252:         C_NATOMS = NATOMS
270:         C_CEIG = CEIG253:         C_CEIG = CEIG
295:         C_BFGSTSTOL = BFGSTSTOL278:         C_BFGSTSTOL = BFGSTSTOL
296:         C_MUPDATE = MUPDATE279:         C_MUPDATE = MUPDATE
297:         C_GMAX = GMAX280:         C_GMAX = GMAX
298:         C_MAXBFGS = MAXBFGS281:         C_MAXBFGS = MAXBFGS
299:         C_MAXERISE = MAXERISE282:         C_MAXERISE = MAXERISE
300:         C_DGUESS = DGUESS283:         C_DGUESS = DGUESS
301:         C_CONVU = CONVU284:         C_CONVU = CONVU
302:         C_FREEZE = FREEZE285:         C_FREEZE = FREEZE
303:         C_NFREEZE = NFREEZE286:         C_NFREEZE = NFREEZE
304:         C_AACONVERGENCET = AACONVERGENCET287:         C_AACONVERGENCET = AACONVERGENCET
305:         C_NADDTARGET = NADDTARGET 
306: 288: 
307:         ! 'C_' prefix denotes those variables which have intent out or inout or are copies of those from common blocks/modules289:         ! 'C_' prefix denotes those variables which have intent out or inout or are copies of those from common blocks/modules
308:         CALL CUDA_BFGSTS(C_NATOMS, C_COORDS, C_CEIG, C_MFLAG, C_ENERGY, C_NEVS, ITMAX, C_ITER, C_MAXXBFGS, C_XMAXERISE, C_RMS, &290:         CALL CUDA_BFGSTS(C_NATOMS, C_COORDS, C_CEIG, C_MFLAG, C_ENERGY, C_NEVS, ITMAX, C_ITER, C_MAXXBFGS, C_XMAXERISE, C_RMS, &
309:                         C_CUDAPOT, C_DEBUG, C_CUDATIMET, NCALLS, C_CFUSIONT, C_COLDFUSIONLIMIT, C_XDGUESS, C_XMUPDATE, C_EVALMIN, &291:                         C_CUDAPOT, C_DEBUG, C_CUDATIMET, NCALLS, C_CFUSIONT, C_COLDFUSIONLIMIT, C_XDGUESS, C_XMUPDATE, C_EVALMIN, &
310:                         C_VECS, C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, &292:                         C_VECS, C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, &
311:                         C_SITESRIGIDBODY, C_RIGIDSINGLES, C_IINVERSE, C_NSECDIAG, C_EVPC, C_PUSHCUT, C_PUSHOFF, &293:                         C_SITESRIGIDBODY, C_RIGIDSINGLES, C_IINVERSE, C_NSECDIAG, C_EVPC, C_PUSHCUT, C_PUSHOFF, &
312:                         C_STPMAX, C_MAXMAX, C_MINMAX, C_TRAD, C_NBFGSMAX1, C_NBFGSMAX2, C_BFGSTSTOL, C_MUPDATE, C_GMAX, &294:                         C_STPMAX, C_MAXMAX, C_MINMAX, C_TRAD, C_NBFGSMAX1, C_NBFGSMAX2, C_BFGSTSTOL, C_MUPDATE, C_GMAX, &
313:                         C_MAXBFGS, C_MAXERISE, C_DGUESS, C_CONVU, C_FREEZE, C_FROZEN, C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET, &295:                         C_MAXBFGS, C_MAXERISE, C_DGUESS, C_CONVU, C_FREEZE, C_FROZEN, C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET)
314:                         C_NADDTARGET, C_LJADDREP, C_LJADDATT) 
315: 296: 
316:         ! Make sure C types with intent out or inout are coverted back to Fortran ones297:         ! Make sure C types with intent out or inout are coverted back to Fortran ones
317: 298: 
318:         DO I = 1,3*NATOMS299:         DO I = 1,3*NATOMS
319:             COORDS(I) = DBLE(C_COORDS(I))300:             COORDS(I) = DBLE(C_COORDS(I))
320:         END DO301:         END DO
321: 302: 
322:         DO I = 1,3*NATOMS303:         DO I = 1,3*NATOMS
323:             VECS(I) = DBLE(C_VECS(I))304:             VECS(I) = DBLE(C_VECS(I))
324:         END DO305:         END DO


r32161/modcudalbfgs.f90 2017-03-22 12:30:31.076280228 +0000 r32160/modcudalbfgs.f90 2017-03-22 12:30:34.632327140 +0000
  1: MODULE MODCUDALBFGS  1: MODULE MODCUDALBFGS
  2:   2: 
  3: USE KEY, ONLY : GMAX, CUDATIMET, CUDAPOT, MAXBFGS, MAXERISE, CFUSIONT, COLDFUSIONLIMIT, DGUESS, GRADSQ, FREEZE, REVERSEUPHILLT, &   3: USE KEY, ONLY : GMAX, CUDATIMET, CUDAPOT, MAXBFGS, MAXERISE, CFUSIONT, COLDFUSIONLIMIT, DGUESS, GRADSQ, FREEZE, REVERSEUPHILLT, & 
  4:                 FREEZE, FROZEN, NFREEZE, PRINTPTS, TWOENDS, PV, DUMPMAG, NSTEPMIN, FIXAFTER, AMBER12T, LJADD3T, NADDTARGET, &   4:                 FREEZE, FROZEN, NFREEZE, PRINTPTS, TWOENDS, PV, DUMPMAG, NSTEPMIN, FIXAFTER, AMBER12T
  5:                 LJADDREP, LJADDATT 
  6:   5: 
  7: USE COMMONS, ONLY : DEBUG  6: USE COMMONS, ONLY : DEBUG
  8:   7: 
  9: USE GENRIGID, ONLY : ATOMRIGIDCOORDT, DEGFREEDOMS, NRIGIDBODY, NSITEPERBODY, RIGIDGROUPS, MAXSITE, SITESRIGIDBODY, &   8: USE GENRIGID, ONLY : ATOMRIGIDCOORDT, DEGFREEDOMS, NRIGIDBODY, NSITEPERBODY, RIGIDGROUPS, MAXSITE, SITESRIGIDBODY, & 
 10:                      RIGIDSINGLES, IINVERSE, AACONVERGENCET  9:                      RIGIDSINGLES, IINVERSE, AACONVERGENCET
 11:  10: 
 12: USE INTCOMMONS, ONLY : BONDSFROMFILE 11: USE INTCOMMONS, ONLY : BONDSFROMFILE
 13:  12: 
 14: USE INTERNALS_WRAPPER, ONLY : INTWRAP_USEINTERNALS 13: USE INTERNALS_WRAPPER, ONLY : INTWRAP_USEINTERNALS
 15:  14: 
 16: USE MODAMBER9, ONLY : CHECKCISTRANSALWAYS, CHECKCISTRANSALWAYSDNA, CHECKCISTRANSALWAYSRNA 15: USE MODAMBER9, ONLY : CHECKCISTRANSALWAYS, CHECKCISTRANSALWAYSDNA, CHECKCISTRANSALWAYSRNA
 17:  16: 
 18: USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_PTR, C_INT, C_DOUBLE, C_BOOL, C_CHAR 17: USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_PTR, C_INT, C_DOUBLE, C_BOOL, C_CHAR
 19:  18: 
 20: IMPLICIT NONE 19: IMPLICIT NONE
 21:  20: 
 22: INTERFACE 21: INTERFACE
 23:     SUBROUTINE CUDA_LBFGS(N, C_XCOORDS, C_GMAX, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, &  22:     SUBROUTINE CUDA_LBFGS(N, C_XCOORDS, C_GMAX, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, & 
 24:                          C_DEBUG, C_CUDATIMET, NCALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, C_DGUESS, MUPDATE, C_ATOMRIGIDCOORDT, &  23:                          C_DEBUG, C_CUDATIMET, NCALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, C_DGUESS, MUPDATE, C_ATOMRIGIDCOORDT, & 
 25:                          C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, C_SITESRIGIDBODY, C_RIGIDSINGLES, &  24:                          C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, C_SITESRIGIDBODY, C_RIGIDSINGLES, & 
 26:                          C_BQMAX, C_IINVERSE, PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET, &  25:                          C_BQMAX, C_IINVERSE, PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET) & 
 27:                          C_NADDTARGET, C_LJADDREP, C_LJADDATT) BIND(C,NAME="setup_lbfgs") 26:                          BIND(C,NAME="setup_lbfgs")
 28:  27: 
 29:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR 28:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR
 30:  29: 
 31:         INTEGER(KIND=C_INT), INTENT(IN) :: N, & ! 3*no. of atoms 30:         INTEGER(KIND=C_INT), INTENT(IN) :: N, & ! 3*no. of atoms
 32:                                            ITMAX, & ! Max. no. of steps allowed in minmization 31:                                            ITMAX, & ! Max. no. of steps allowed in minmization
 33:                                            C_DEGFREEDOMS, & ! Rigid Body Framework (RBF): no. of degrees of freedom 32:                                            C_DEGFREEDOMS, & ! Rigid Body Framework (RBF): no. of degrees of freedom
 34:                                            C_NRIGIDBODY, & ! RBF: no. of rigid bodies 33:                                            C_NRIGIDBODY, & ! RBF: no. of rigid bodies
 35:                                            C_MAXSITE, & ! RBF: max. no. of sites in a rigid body 34:                                            C_MAXSITE, & ! RBF: max. no. of sites in a rigid body
 36:                                            MUPDATE, & ! History size 35:                                            MUPDATE, & ! History size
 37:                                            C_NFREEZE, & ! No. of frozen atoms 36:                                            C_NFREEZE ! No. of frozen atoms
 38:                                            C_NADDTARGET ! Target cluster size (addressability) 
 39:  37: 
 40:         INTEGER(KIND=C_INT), DIMENSION(C_NRIGIDBODY), INTENT(IN) :: C_NSITEPERBODY ! RBF: no. of rigid body sites 38:         INTEGER(KIND=C_INT), DIMENSION(C_NRIGIDBODY), INTENT(IN) :: C_NSITEPERBODY ! RBF: no. of rigid body sites
 41:  39: 
 42:         INTEGER(KIND=C_INT), DIMENSION(C_MAXSITE*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDGROUPS ! RBF: list of atoms in rigid bodies 40:         INTEGER(KIND=C_INT), DIMENSION(C_MAXSITE*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDGROUPS ! RBF: list of atoms in rigid bodies
 43:  41: 
 44:  42: 
 45:         INTEGER(KIND=C_INT), DIMENSION(C_DEGFREEDOMS/3 - 2*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDSINGLES ! RBF: list of atoms not in rigid bodies 43:         INTEGER(KIND=C_INT), DIMENSION(C_DEGFREEDOMS/3 - 2*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDSINGLES ! RBF: list of atoms not in rigid bodies
 46:  44: 
 47:         INTEGER(KIND=C_INT), INTENT(OUT) :: C_ITDONE, & ! No. of LBFGS iterations done 45:         INTEGER(KIND=C_INT), INTENT(OUT) :: C_ITDONE, & ! No. of LBFGS iterations done
 48:                                             NCALLS ! No. of potential calls made during this call to LBFGS 46:                                             NCALLS ! No. of potential calls made during this call to LBFGS
 51:                                            C_MAXBFGS, & ! Max. step size allowed  49:                                            C_MAXBFGS, & ! Max. step size allowed 
 52:                                            C_MAXERISE, & ! Max. energy rise allowed  50:                                            C_MAXERISE, & ! Max. energy rise allowed 
 53:                                            C_COLDFUSIONLIMIT, & ! Limit below which cold fusion is diagnosed and minimization terminated 51:                                            C_COLDFUSIONLIMIT, & ! Limit below which cold fusion is diagnosed and minimization terminated
 54:                                            C_DGUESS, & ! Initial guess for inverse Hessian diagonal elements 52:                                            C_DGUESS, & ! Initial guess for inverse Hessian diagonal elements
 55:                                            C_BQMAX ! (1/5) of threshold under which aaconvergence is used, equivalent to C_GMAX 53:                                            C_BQMAX ! (1/5) of threshold under which aaconvergence is used, equivalent to C_GMAX
 56:  54: 
 57:         REAL(KIND=C_DOUBLE), DIMENSION(C_MAXSITE*3*C_NRIGIDBODY), INTENT(IN) :: C_SITESRIGIDBODY ! RBF: coordinates of the rigid body sites 55:         REAL(KIND=C_DOUBLE), DIMENSION(C_MAXSITE*3*C_NRIGIDBODY), INTENT(IN) :: C_SITESRIGIDBODY ! RBF: coordinates of the rigid body sites
 58:  56: 
 59:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration 57:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration
 60:  58: 
 61:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 
 62:  
 63:         REAL(KIND=C_DOUBLE), DIMENSION(N), INTENT(INOUT) :: C_XCOORDS ! Coordinates 59:         REAL(KIND=C_DOUBLE), DIMENSION(N), INTENT(INOUT) :: C_XCOORDS ! Coordinates
 64:  60: 
 65:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy 61:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy
 66:                                             C_RMS, & ! RMS force 62:                                             C_RMS, & ! RMS force
 67:                                             POTENTIALTIME ! Time taken in calculating potential 63:                                             POTENTIALTIME ! Time taken in calculating potential
 68:  64: 
 69:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info. 65:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info.
 70:                                             C_CUDATIMET, & ! If true, print timing info. 66:                                             C_CUDATIMET, & ! If true, print timing info.
 71:                                             C_ATOMRIGIDCOORDT, & ! If false, use rigid body coordinates 67:                                             C_ATOMRIGIDCOORDT, & ! If false, use rigid body coordinates
 72:                                             PROJECT, & ! If true, project out the components of the gradient along the eigenvector 68:                                             PROJECT, & ! If true, project out the components of the gradient along the eigenvector
 77:  73: 
 78:         LOGICAL(KIND=C_BOOL), INTENT(OUT) :: C_MFLAG, & ! True if quench converged 74:         LOGICAL(KIND=C_BOOL), INTENT(OUT) :: C_MFLAG, & ! True if quench converged
 79:                                              C_COLDFUSION ! Set to true during minimization if cold fusion diagnosed 75:                                              C_COLDFUSION ! Set to true during minimization if cold fusion diagnosed
 80:  76: 
 81:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used 77:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used
 82:  78: 
 83:     END SUBROUTINE CUDA_LBFGS 79:     END SUBROUTINE CUDA_LBFGS
 84: END INTERFACE 80: END INTERFACE
 85:  81: 
 86: INTERFACE 82: INTERFACE
 87:     SUBROUTINE CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, &  83:     SUBROUTINE CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS) BIND(C,NAME="gminoptim_enegrad_cputogpu")
 88:                                      C_LJADDATT) BIND(C,NAME="gminoptim_enegrad_cputogpu") 
 89:  84: 
 90:         IMPORT :: C_INT, C_DOUBLE 85:         IMPORT :: C_INT, C_DOUBLE
 91:  86: 
 92:         INTEGER(KIND=C_INT), INTENT(IN) :: NATOMS, & ! No. of atoms 87:         INTEGER(KIND=C_INT), INTENT(IN) :: NATOMS ! No. of atoms
 93:                                            C_NADDTARGET ! Target cluster size (addressability) 
 94:  88: 
 95:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(IN) :: COORDS ! Atomic coordinates 89:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(IN) :: COORDS ! Atomic coordinates
 96:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 
 97:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_TOTENERGY ! Total energy of the system 90:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_TOTENERGY ! Total energy of the system
 98:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(OUT) :: C_GRADIENTS ! Gradient of the energy w.r.t. each atomic coordinate 91:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(OUT) :: C_GRADIENTS ! Gradient of the energy w.r.t. each atomic coordinate
 99:  92: 
100:     END SUBROUTINE CUDA_ENEGRAD_CPUTOGPU 93:     END SUBROUTINE CUDA_ENEGRAD_CPUTOGPU
101: END INTERFACE 94: END INTERFACE
102:  95: 
103: CONTAINS 96: CONTAINS
104:  97: 
105:     SUBROUTINE CUDA_LBFGS_WRAPPER(N, MUPDATE, XCOORDS, MFLAG, ENERGY, RMS, ITMAX, ITDONE, PROJECT) 98:     SUBROUTINE CUDA_LBFGS_WRAPPER(N, MUPDATE, XCOORDS, MFLAG, ENERGY, RMS, ITMAX, ITDONE, PROJECT)
106:  99: 
107:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_LBFGS are converted directly100:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_LBFGS are converted directly
108:         INTEGER(KIND=C_INT) :: N, ITMAX, C_ITDONE, NCALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, MUPDATE, C_NFREEZE, C_NADDTARGET101:         INTEGER(KIND=C_INT) :: N, ITMAX, C_ITDONE, NCALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, MUPDATE, C_NFREEZE
109:         INTEGER(KIND=C_INT), DIMENSION(NRIGIDBODY) :: C_NSITEPERBODY102:         INTEGER(KIND=C_INT), DIMENSION(NRIGIDBODY) :: C_NSITEPERBODY
110:         INTEGER(KIND=C_INT), DIMENSION(MAXSITE*NRIGIDBODY) :: C_RIGIDGROUPS103:         INTEGER(KIND=C_INT), DIMENSION(MAXSITE*NRIGIDBODY) :: C_RIGIDGROUPS
111:         INTEGER(KIND=C_INT), DIMENSION(DEGFREEDOMS/3 - 2*NRIGIDBODY) :: C_RIGIDSINGLES104:         INTEGER(KIND=C_INT), DIMENSION(DEGFREEDOMS/3 - 2*NRIGIDBODY) :: C_RIGIDSINGLES
112: 105: 
113:         REAL(KIND=C_DOUBLE) :: C_MAXBFGS, C_MAXERISE, C_ENERGY, C_RMS, C_COLDFUSIONLIMIT, C_DGUESS, C_BQMAX, C_GMAX, POTENTIALTIME106:         REAL(KIND=C_DOUBLE) :: C_MAXBFGS, C_MAXERISE, C_ENERGY, C_RMS, C_COLDFUSIONLIMIT, C_DGUESS, C_BQMAX, C_GMAX, POTENTIALTIME
114:         REAL(KIND=C_DOUBLE), DIMENSION(MAXSITE*3*NRIGIDBODY) :: C_SITESRIGIDBODY107:         REAL(KIND=C_DOUBLE), DIMENSION(MAXSITE*3*NRIGIDBODY) :: C_SITESRIGIDBODY
115:         REAL(KIND=C_DOUBLE), DIMENSION(NRIGIDBODY*3*3) :: C_IINVERSE108:         REAL(KIND=C_DOUBLE), DIMENSION(NRIGIDBODY*3*3) :: C_IINVERSE
116:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT 
117:         REAL(KIND=C_DOUBLE), DIMENSION(N) :: C_XCOORDS109:         REAL(KIND=C_DOUBLE), DIMENSION(N) :: C_XCOORDS
118: 110: 
119:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_ATOMRIGIDCOORDT, C_COLDFUSION, PROJECT, C_FREEZE, C_AACONVERGENCET111:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_ATOMRIGIDCOORDT, C_COLDFUSION, PROJECT, C_FREEZE, C_AACONVERGENCET
120:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3) :: C_FROZEN112:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3) :: C_FROZEN
121: 113: 
122:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT114:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
123: 115: 
124:         ! Variables passed as *arguments through this wrapper* (not common) with intent out for CUDA_LBFGS are not passed into it116:         ! Variables passed as *arguments through this wrapper* (not common) with intent out for CUDA_LBFGS are not passed into it
125:         ! Therefore uninitialised C types are passed in and converted types are copied back after the call117:         ! Therefore uninitialised C types are passed in and converted types are copied back after the call
126: 118: 
182:         END DO174:         END DO
183: 175: 
184:         DO K = 1,3176:         DO K = 1,3
185:             DO J = 1,3177:             DO J = 1,3
186:                 DO I = 1,NRIGIDBODY178:                 DO I = 1,NRIGIDBODY
187:                    C_IINVERSE((K - 1)*3*NRIGIDBODY + (J - 1)*NRIGIDBODY + I) = IINVERSE(I,J,K)179:                    C_IINVERSE((K - 1)*3*NRIGIDBODY + (J - 1)*NRIGIDBODY + I) = IINVERSE(I,J,K)
188:                 END DO180:                 END DO
189:             END DO181:             END DO
190:         END DO182:         END DO
191: 183: 
192:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN 
193:             DO J = 1,NADDTARGET 
194:                 DO I = 1,NADDTARGET 
195:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J) 
196:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J) 
197:                 END DO 
198:             END DO 
199:         ELSE 
200:             C_LJADDREP(:) = 1.0 
201:             C_LJADDATT(:) = 1.0 
202:         END IF 
203:  
204:         DO I = 1,N184:         DO I = 1,N
205:             C_XCOORDS(I) = XCOORDS(I)185:             C_XCOORDS(I) = XCOORDS(I)
206:         END DO186:         END DO
207: 187: 
208:         DO I = 1,(N/3)188:         DO I = 1,(N/3)
209:             C_FROZEN(I) = FROZEN(I)189:             C_FROZEN(I) = FROZEN(I)
210:         END DO190:         END DO
211: 191: 
212:         C_CUDAPOT = CUDAPOT192:         C_CUDAPOT = CUDAPOT
213:         C_DEBUG = DEBUG193:         C_DEBUG = DEBUG
218:         C_DEGFREEDOMS = DEGFREEDOMS198:         C_DEGFREEDOMS = DEGFREEDOMS
219:         C_NRIGIDBODY = NRIGIDBODY199:         C_NRIGIDBODY = NRIGIDBODY
220:         C_MAXSITE = MAXSITE200:         C_MAXSITE = MAXSITE
221:         C_COLDFUSIONLIMIT = COLDFUSIONLIMIT201:         C_COLDFUSIONLIMIT = COLDFUSIONLIMIT
222:         C_DGUESS = DGUESS202:         C_DGUESS = DGUESS
223:         C_BQMAX = GMAX203:         C_BQMAX = GMAX
224:         C_GMAX = GMAX204:         C_GMAX = GMAX
225:         C_FREEZE = FREEZE205:         C_FREEZE = FREEZE
226:         C_NFREEZE = NFREEZE206:         C_NFREEZE = NFREEZE
227:         C_AACONVERGENCET = AACONVERGENCET207:         C_AACONVERGENCET = AACONVERGENCET
228:         C_NADDTARGET = NADDTARGET 
229: 208: 
230:         ! 'C_' prefix denotes those variables which have intent out or inout or are copies of those from common blocks/modules 209:         ! 'C_' prefix denotes those variables which have intent out or inout or are copies of those from common blocks/modules 
231:         CALL CUDA_LBFGS(N, C_XCOORDS, C_GMAX, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, & 210:         CALL CUDA_LBFGS(N, C_XCOORDS, C_GMAX, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, & 
232:                        C_DEBUG, C_CUDATIMET, NCALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, C_DGUESS, MUPDATE, & 211:                        C_DEBUG, C_CUDATIMET, NCALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, C_DGUESS, MUPDATE, & 
233:                        C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, & 212:                        C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, & 
234:                        C_SITESRIGIDBODY, C_RIGIDSINGLES, C_BQMAX, C_IINVERSE, PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, & 213:                        C_SITESRIGIDBODY, C_RIGIDSINGLES, C_BQMAX, C_IINVERSE, PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, & 
235:                        POTENTIALTIME, C_AACONVERGENCET, C_NADDTARGET, C_LJADDREP, C_LJADDATT)214:                        POTENTIALTIME, C_AACONVERGENCET)
236: 215: 
237:         ! Make sure C types with intent out or inout are coverted back to Fortran ones216:         ! Make sure C types with intent out or inout are coverted back to Fortran ones
238: 217: 
239:         DO I = 1,N218:         DO I = 1,N
240:             XCOORDS(I) = DBLE(C_XCOORDS(I))219:             XCOORDS(I) = DBLE(C_XCOORDS(I))
241:         END DO220:         END DO
242: 221: 
243:         ENERGY = DBLE(C_ENERGY)222:         ENERGY = DBLE(C_ENERGY)
244:         RMS = DBLE(C_RMS)223:         RMS = DBLE(C_RMS)
245:         MFLAG = LOGICAL(C_MFLAG)224:         MFLAG = LOGICAL(C_MFLAG)
257: 236: 
258:         KNOWE=.FALSE.237:         KNOWE=.FALSE.
259:         KNOWG=.FALSE.238:         KNOWG=.FALSE.
260:         KNOWH=.FALSE.239:         KNOWH=.FALSE.
261: 240: 
262:     END SUBROUTINE CUDA_LBFGS_WRAPPER241:     END SUBROUTINE CUDA_LBFGS_WRAPPER
263: 242: 
264: 243: 
265:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)244:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)
266: 245: 
267:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET246:         INTEGER(KIND=C_INT) :: NATOMS
268: 247: 
269:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY248:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY
270:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS249:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS
271:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT 
272: 250: 
273:         INTEGER :: X, I, J251:         INTEGER :: X
274: 252: 
275:         DOUBLE PRECISION :: TOTENERGY253:         DOUBLE PRECISION :: TOTENERGY
276:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS254:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS
277: 255: 
278:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN 
279:             DO J = 1,NADDTARGET 
280:                 DO I = 1,NADDTARGET 
281:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J) 
282:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J) 
283:                 END DO 
284:             END DO 
285:         ELSE 
286:             C_LJADDREP(:) = 1.0 
287:             C_LJADDATT(:) = 1.0 
288:         END IF 
289:  
290:         C_NADDTARGET = NADDTARGET 
291:  
292:         ! Calculates the energy and gradients on the GPU using the GB potential256:         ! Calculates the energy and gradients on the GPU using the GB potential
293:         CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)257:         CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS)
294: 258: 
295:         TOTENERGY = DBLE(C_TOTENERGY)259:         TOTENERGY = DBLE(C_TOTENERGY)
296: 260: 
297:         DO X = 1,(3*NATOMS)261:         DO X = 1,(3*NATOMS)
298:             GRADIENTS(X) = DBLE(C_GRADIENTS(X))262:             GRADIENTS(X) = DBLE(C_GRADIENTS(X))
299:         END DO263:         END DO
300: 264: 
301:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER265:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER
302: 266: 
303: END MODULE MODCUDALBFGS267: END MODULE MODCUDALBFGS


r32161/modcudalbfgs.F90 2017-03-22 12:30:29.748262703 +0000 r32160/modcudalbfgs.F90 2017-03-22 12:30:33.292309460 +0000
  1: MODULE MODCUDALBFGS  1: MODULE MODCUDALBFGS
  2:   2: 
  3: #ifndef DUMMY_CUDA  3: #ifndef DUMMY_CUDA
  4:   4: 
  5: USE COMMONS, ONLY : RMS, DEBUG, CUDATIMET, MAXBFGS, MAXERISE, CUDAPOT, NPCALL, COLDFUSION, COLDFUSIONLIMIT, MYUNIT, DGUESS, &  5: USE COMMONS, ONLY : RMS, DEBUG, CUDATIMET, MAXBFGS, MAXERISE, CUDAPOT, NPCALL, COLDFUSION, COLDFUSIONLIMIT, MYUNIT, DGUESS, &
  6:                     BQMAX, FREEZE, SEEDT, FREEZECORE, FROZEN, NFREEZE, QUENCHDOS, CENT, CALCQT, INTMINT, DUMPT, COMPRESSRIGIDT, &   6:                     BQMAX, FREEZE, SEEDT, FREEZECORE, FROZEN, NFREEZE, QUENCHDOS, CENT, CALCQT, INTMINT, DUMPT, COMPRESSRIGIDT, & 
  7:                     AMBER12T, LJADD3T, NADDTARGET, LJADDREP, LJADDATT  7:                     AMBER12T
  8:   8: 
  9: USE GENRIGID, ONLY : ATOMRIGIDCOORDT, DEGFREEDOMS, NRIGIDBODY, NSITEPERBODY, RIGIDGROUPS, MAXSITE, SITESRIGIDBODY, &   9: USE GENRIGID, ONLY : ATOMRIGIDCOORDT, DEGFREEDOMS, NRIGIDBODY, NSITEPERBODY, RIGIDGROUPS, MAXSITE, SITESRIGIDBODY, & 
 10:                      RIGIDSINGLES, IINVERSE, AACONVERGENCET 10:                      RIGIDSINGLES, IINVERSE, AACONVERGENCET
 11:  11: 
 12: USE MODAMBER9, ONLY : STEEREDMINT 12: USE MODAMBER9, ONLY : STEEREDMINT
 13:  13: 
 14: #endif 14: #endif
 15:  15: 
 16: USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_PTR, C_INT, C_DOUBLE, C_BOOL, C_CHAR 16: USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_PTR, C_INT, C_DOUBLE, C_BOOL, C_CHAR
 17:  17: 
 18: IMPLICIT NONE 18: IMPLICIT NONE
 19:  19: 
 20: #ifndef DUMMY_CUDA 20: #ifndef DUMMY_CUDA
 21: INTERFACE 21: INTERFACE
 22:     SUBROUTINE CUDA_LBFGS(N, C_XCOORDS, EPS, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, &  22:     SUBROUTINE CUDA_LBFGS(N, C_XCOORDS, EPS, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, & 
 23:                          C_DEBUG, C_CUDATIMET, ECALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, & 23:                          C_DEBUG, C_CUDATIMET, ECALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, &
 24:                          C_DGUESS, MUPDATE, C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, & 24:                          C_DGUESS, MUPDATE, C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, &
 25:                          C_RIGIDGROUPS, C_MAXSITE, C_SITESRIGIDBODY, C_RIGIDSINGLES, C_BQMAX, C_IINVERSE, & 25:                          C_RIGIDGROUPS, C_MAXSITE, C_SITESRIGIDBODY, C_RIGIDSINGLES, C_BQMAX, C_IINVERSE, &
 26:                          PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET, C_NADDTARGET, C_LJADDREP, &  26:                          PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, POTENTIALTIME, C_AACONVERGENCET) BIND(C,NAME="setup_lbfgs")
 27:                          C_LJADDATT) BIND(C,NAME="setup_lbfgs") 
 28:  27: 
 29:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR 28:         IMPORT :: C_INT, C_DOUBLE, C_BOOL, C_CHAR
 30:  29: 
 31:         INTEGER(KIND=C_INT), INTENT(IN) :: N, & ! 3*no. of atoms 30:         INTEGER(KIND=C_INT), INTENT(IN) :: N, & ! 3*no. of atoms
 32:                                            ITMAX, & ! Max. no. of steps allowed in minmization 31:                                            ITMAX, & ! Max. no. of steps allowed in minmization
 33:                                            C_DEGFREEDOMS, & ! Rigid Body Framework (RBF): no. of degrees of freedom 32:                                            C_DEGFREEDOMS, & ! Rigid Body Framework (RBF): no. of degrees of freedom
 34:                                            C_NRIGIDBODY, & ! RBF: no. of rigid bodies 33:                                            C_NRIGIDBODY, & ! RBF: no. of rigid bodies
 35:                                            C_MAXSITE, & ! RBF: max. no. of sites in a rigid body 34:                                            C_MAXSITE, & ! RBF: max. no. of sites in a rigid body
 36:                                            MUPDATE, & ! History size 35:                                            MUPDATE, & ! History size
 37:                                            C_NFREEZE, & ! No. of frozen atoms 36:                                            C_NFREEZE ! No. of frozen atoms
 38:                                            C_NADDTARGET ! Target cluster size (addressability) 
 39:  37: 
 40:         INTEGER(KIND=C_INT), DIMENSION(C_NRIGIDBODY), INTENT(IN) :: C_NSITEPERBODY ! RBF: no. of rigid body sites 38:         INTEGER(KIND=C_INT), DIMENSION(C_NRIGIDBODY), INTENT(IN) :: C_NSITEPERBODY ! RBF: no. of rigid body sites
 41:  39: 
 42:         INTEGER(KIND=C_INT), DIMENSION(C_MAXSITE*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDGROUPS ! RBF: list of atoms in rigid bodies 40:         INTEGER(KIND=C_INT), DIMENSION(C_MAXSITE*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDGROUPS ! RBF: list of atoms in rigid bodies
 43:  41: 
 44:         INTEGER(KIND=C_INT), DIMENSION(C_DEGFREEDOMS/3 - 2*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDSINGLES ! RBF: list of atoms not in rigid bodies 42:         INTEGER(KIND=C_INT), DIMENSION(C_DEGFREEDOMS/3 - 2*C_NRIGIDBODY), INTENT(IN) :: C_RIGIDSINGLES ! RBF: list of atoms not in rigid bodies
 45:  43: 
 46:         INTEGER(KIND=C_INT), INTENT(OUT) :: C_ITDONE, & ! No. of LBFGS iterations done 44:         INTEGER(KIND=C_INT), INTENT(OUT) :: C_ITDONE, & ! No. of LBFGS iterations done
 47:                                             ECALLS ! Number of potential calls made during this call to LBFGS 45:                                             ECALLS ! Number of potential calls made during this call to LBFGS
 48:  46: 
 50:                                            C_MAXBFGS, & ! Max. step size allowed 48:                                            C_MAXBFGS, & ! Max. step size allowed
 51:                                            C_MAXERISE, & ! Max. energy rise allowed 49:                                            C_MAXERISE, & ! Max. energy rise allowed
 52:                                            C_COLDFUSIONLIMIT, & ! Limit below which cold fusion is diagnosed and minimization terminated 50:                                            C_COLDFUSIONLIMIT, & ! Limit below which cold fusion is diagnosed and minimization terminated
 53:                                            C_DGUESS, & ! Initial guess for inverse Hessian diagonal elements 51:                                            C_DGUESS, & ! Initial guess for inverse Hessian diagonal elements
 54:                                            C_BQMAX ! Sloppy quench tolerance for RMS gradient 52:                                            C_BQMAX ! Sloppy quench tolerance for RMS gradient
 55:  53: 
 56:         REAL(KIND=C_DOUBLE), DIMENSION(C_MAXSITE*3*C_NRIGIDBODY), INTENT(IN) :: C_SITESRIGIDBODY ! RBF: coordinates of the rigid body sites 54:         REAL(KIND=C_DOUBLE), DIMENSION(C_MAXSITE*3*C_NRIGIDBODY), INTENT(IN) :: C_SITESRIGIDBODY ! RBF: coordinates of the rigid body sites
 57:  55: 
 58:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration 56:         REAL(KIND=C_DOUBLE), DIMENSION(C_NRIGIDBODY*3*3), INTENT(IN) :: C_IINVERSE ! RBF: inverse eigenvalues of the unweighted tensor of gyration
 59:  57: 
 60:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 
 61:  
 62:         REAL(KIND=C_DOUBLE), DIMENSION(N), INTENT(INOUT) :: C_XCOORDS ! Coordinates 58:         REAL(KIND=C_DOUBLE), DIMENSION(N), INTENT(INOUT) :: C_XCOORDS ! Coordinates
 63:  59: 
 64:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy 60:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_ENERGY, & ! Energy
 65:                                             C_RMS, & ! RMS force 61:                                             C_RMS, & ! RMS force
 66:                                             POTENTIALTIME ! Time taken in calculating potential - not used in GMIN 62:                                             POTENTIALTIME ! Time taken in calculating potential - not used in GMIN
 67:  63: 
 68:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info. 64:         LOGICAL(KIND=C_BOOL), INTENT(IN) :: C_DEBUG, & ! If true, print debug info.
 69:                                             C_CUDATIMET, & ! If true, print timing info.  65:                                             C_CUDATIMET, & ! If true, print timing info. 
 70:                                             C_ATOMRIGIDCOORDT, & ! If false, use rigid body coordinates 66:                                             C_ATOMRIGIDCOORDT, & ! If false, use rigid body coordinates
 71:                                             C_FREEZE, & ! If true, freeze some specified atoms 67:                                             C_FREEZE, & ! If true, freeze some specified atoms
 76:  72: 
 77:         LOGICAL(KIND=C_BOOL), INTENT(OUT) :: C_MFLAG, & ! True if quench converged 73:         LOGICAL(KIND=C_BOOL), INTENT(OUT) :: C_MFLAG, & ! True if quench converged
 78:                                              C_COLDFUSION ! Set to true during minimization if cold fusion diagnosed 74:                                              C_COLDFUSION ! Set to true during minimization if cold fusion diagnosed
 79:  75: 
 80:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used 76:         CHARACTER(LEN=1, KIND=C_CHAR), INTENT(IN) :: C_CUDAPOT ! Character specifying the CUDA potential to be used
 81:  77: 
 82:     END SUBROUTINE CUDA_LBFGS 78:     END SUBROUTINE CUDA_LBFGS
 83: END INTERFACE 79: END INTERFACE
 84:  80: 
 85: INTERFACE 81: INTERFACE
 86:     SUBROUTINE CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, &  82:     SUBROUTINE CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS) BIND(C,NAME="gminoptim_enegrad_cputogpu")
 87:                                      C_LJADDATT) BIND(C,NAME="gminoptim_enegrad_cputogpu") 
 88:  83: 
 89:         IMPORT :: C_INT, C_DOUBLE 84:         IMPORT :: C_INT, C_DOUBLE
 90:  85: 
 91:         INTEGER(KIND=C_INT), INTENT(IN) :: NATOMS, & ! No. of atoms 86:         INTEGER(KIND=C_INT), INTENT(IN) :: NATOMS ! No. of atoms
 92:                                            C_NADDTARGET ! Target cluster size (addressability) 
 93:  87: 
 94:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(IN) :: COORDS ! Atomic coordinates 88:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(IN) :: COORDS ! Atomic coordinates
 95:         REAL(KIND=C_DOUBLE), DIMENSION(C_NADDTARGET, C_NADDTARGET), INTENT(IN) ::C_LJADDREP, C_LJADDATT ! Repulsive/attractive epsilon matrix 
 96:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_TOTENERGY ! Total energy of the system 89:         REAL(KIND=C_DOUBLE), INTENT(OUT) :: C_TOTENERGY ! Total energy of the system
 97:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(OUT) :: C_GRADIENTS ! Gradient of the energy w.r.t. each atomic coordinate 90:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS), INTENT(OUT) :: C_GRADIENTS ! Gradient of the energy w.r.t. each atomic coordinate
 98:  91: 
 99:     END SUBROUTINE CUDA_ENEGRAD_CPUTOGPU     92:     END SUBROUTINE CUDA_ENEGRAD_CPUTOGPU    
100: END INTERFACE 93: END INTERFACE
101: #endif /* DUMMY_CUDA */ 94: #endif /* DUMMY_CUDA */
102:  95: 
103: CONTAINS 96: CONTAINS
104:  97: 
105:     SUBROUTINE CUDA_LBFGS_WRAPPER(N, MUPDATE, XCOORDS, EPS, MFLAG, ENERGY, ITMAX, ITDONE, RESET) 98:     SUBROUTINE CUDA_LBFGS_WRAPPER(N, MUPDATE, XCOORDS, EPS, MFLAG, ENERGY, ITMAX, ITDONE, RESET)
106:  99: 
107:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_LBFGS are converted directly100:         ! Variables passed as *arguments through this wrapper* (not common) with intent in for CUDA_LBFGS are converted directly
108: #ifndef DUMMY_CUDA101: #ifndef DUMMY_CUDA
109:         INTEGER(KIND=C_INT) :: N, MUPDATE, ITMAX, C_ITDONE, ECALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, C_NFREEZE, C_NADDTARGET102:         INTEGER(KIND=C_INT) :: N, MUPDATE, ITMAX, C_ITDONE, ECALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, C_NFREEZE
110:         INTEGER(KIND=C_INT), DIMENSION(NRIGIDBODY) :: C_NSITEPERBODY103:         INTEGER(KIND=C_INT), DIMENSION(NRIGIDBODY) :: C_NSITEPERBODY
111:         INTEGER(KIND=C_INT), DIMENSION(MAXSITE*NRIGIDBODY) :: C_RIGIDGROUPS104:         INTEGER(KIND=C_INT), DIMENSION(MAXSITE*NRIGIDBODY) :: C_RIGIDGROUPS
112:         INTEGER(KIND=C_INT), DIMENSION(DEGFREEDOMS/3 - 2*NRIGIDBODY) :: C_RIGIDSINGLES105:         INTEGER(KIND=C_INT), DIMENSION(DEGFREEDOMS/3 - 2*NRIGIDBODY) :: C_RIGIDSINGLES
113: 106: 
114:         REAL(KIND=C_DOUBLE) :: EPS, C_MAXBFGS, C_MAXERISE, C_ENERGY, C_RMS, C_COLDFUSIONLIMIT, C_DGUESS, C_BQMAX, POTENTIALTIME107:         REAL(KIND=C_DOUBLE) :: EPS, C_MAXBFGS, C_MAXERISE, C_ENERGY, C_RMS, C_COLDFUSIONLIMIT, C_DGUESS, C_BQMAX, POTENTIALTIME
115:         REAL(KIND=C_DOUBLE), DIMENSION(MAXSITE*3*NRIGIDBODY) :: C_SITESRIGIDBODY108:         REAL(KIND=C_DOUBLE), DIMENSION(MAXSITE*3*NRIGIDBODY) :: C_SITESRIGIDBODY
116:         REAL(KIND=C_DOUBLE), DIMENSION(NRIGIDBODY*3*3) :: C_IINVERSE109:         REAL(KIND=C_DOUBLE), DIMENSION(NRIGIDBODY*3*3) :: C_IINVERSE
117:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT 
118:         REAL(KIND=C_DOUBLE), DIMENSION(N) :: C_XCOORDS110:         REAL(KIND=C_DOUBLE), DIMENSION(N) :: C_XCOORDS
119: 111: 
120:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_ATOMRIGIDCOORDT, C_COLDFUSION, C_FREEZE, PROJECT, C_AACONVERGENCET112:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_ATOMRIGIDCOORDT, C_COLDFUSION, C_FREEZE, PROJECT, C_AACONVERGENCET
121:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3) :: C_FROZEN113:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3) :: C_FROZEN
122: 114: 
123:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT115:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
124: #endif116: #endif
125:         ! Same as above, but now dimension(1).117:         ! Same as above, but now dimension(1).
126: #ifdef DUMMY_CUDA118: #ifdef DUMMY_CUDA
127:         INTEGER(KIND=C_INT) :: N, MUPDATE, ITMAX, C_ITDONE, ECALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, C_NFREEZE, C_NADDTARGET119:         INTEGER(KIND=C_INT) :: N, MUPDATE, ITMAX, C_ITDONE, ECALLS, C_DEGFREEDOMS, C_NRIGIDBODY, C_MAXSITE, C_NFREEZE
128:         INTEGER(KIND=C_INT), DIMENSION(1) :: C_NSITEPERBODY120:         INTEGER(KIND=C_INT), DIMENSION(1) :: C_NSITEPERBODY
129:         INTEGER(KIND=C_INT), DIMENSION(1) :: C_RIGIDGROUPS121:         INTEGER(KIND=C_INT), DIMENSION(1) :: C_RIGIDGROUPS
130:         INTEGER(KIND=C_INT), DIMENSION(1) :: C_RIGIDSINGLES122:         INTEGER(KIND=C_INT), DIMENSION(1) :: C_RIGIDSINGLES
131: 123: 
132:         REAL(KIND=C_DOUBLE) :: EPS, C_MAXBFGS, C_MAXERISE, C_ENERGY, C_RMS, C_COLDFUSIONLIMIT, C_DGUESS, C_BQMAX, POTENTIALTIME124:         REAL(KIND=C_DOUBLE) :: EPS, C_MAXBFGS, C_MAXERISE, C_ENERGY, C_RMS, C_COLDFUSIONLIMIT, C_DGUESS, C_BQMAX, POTENTIALTIME
133:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_SITESRIGIDBODY125:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_SITESRIGIDBODY
134:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_IINVERSE126:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_IINVERSE
135:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_LJADDREP, C_LJADDATT 
136:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_XCOORDS127:         REAL(KIND=C_DOUBLE), DIMENSION(1) :: C_XCOORDS
137: 128: 
138:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_ATOMRIGIDCOORDT, C_COLDFUSION, C_FREEZE, PROJECT, C_AACONVERGENCET129:         LOGICAL(KIND=C_BOOL) :: C_MFLAG, C_DEBUG, C_CUDATIMET, C_ATOMRIGIDCOORDT, C_COLDFUSION, C_FREEZE, PROJECT, C_AACONVERGENCET
139:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3) :: C_FROZEN130:         LOGICAL(KIND=C_BOOL), DIMENSION(N/3) :: C_FROZEN
140: 131: 
141:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT132:         CHARACTER(LEN=1, KIND=C_CHAR) :: C_CUDAPOT
142: #endif133: #endif
143: 134: 
144:         ! Variables passed as *arguments through this wrapper* (not common) with intent out for CUDA_LBFGS are not passed into it135:         ! Variables passed as *arguments through this wrapper* (not common) with intent out for CUDA_LBFGS are not passed into it
145:         ! Therefore uninitialised C types are passed in and converted types are copied back after the call136:         ! Therefore uninitialised C types are passed in and converted types are copied back after the call
152:         COMMON /MYPOT/ POTEL143:         COMMON /MYPOT/ POTEL
153: 144: 
154:         IF (.NOT. RESET) THEN145:         IF (.NOT. RESET) THEN
155:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Warning: LBFGS resetting, though RESET is false. "146:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Warning: LBFGS resetting, though RESET is false. "
156:         END IF147:         END IF
157: 148: 
158:         IF (DUMPT .AND. DEBUG) THEN149:         IF (DUMPT .AND. DEBUG) THEN
159:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Warning: printing behaviour of DUMP and DEBUG during minimization is not implemented. "150:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Warning: printing behaviour of DUMP and DEBUG during minimization is not implemented. "
160:         END IF151:         END IF
161: 152: 
162:         IF ((SEEDT.AND.FREEZECORE) .OR. STEEREDMINT .OR. QUENCHDOS .OR. CALCQT .OR. INTMINT .OR. COMPRESSRIGIDT) THEN153:         IF ((SEEDT.AND.FREEZECORE) .OR. STEEREDMINT .OR. QUENCHDOS .OR. CENT .OR. CALCQT .OR. INTMINT .OR. COMPRESSRIGIDT) THEN
163:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Keyword SEED with FREEZECORE is not yet supported. SEED can be used with NOFREEZE. &154:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Keyword SEED with FREEZECORE is not yet supported. SEED can be used with NOFREEZE. &
164:                                  Keywords STEEREDMIN, QUENCHDOS, CALCQ, INTMIN, COMPRESSRIGID are not yet supported. & 155:                                  Keywords STEEREDMIN, QUENCHDOS, CENTRE, CALCQ, INTMIN, COMPRESSRIGID are not yet supported. & 
165:                                  Contact rgm38 if you would like a feature to be added. "156:                                  Contact rgm38 if you would like a feature to be added. "
166:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Disclaimer: this list might not be exhaustive! "157:             WRITE(MYUNIT,'(A)') "modcudalbfgs> Disclaimer: this list might not be exhaustive! "
167:             STOP158:             STOP
168:         END IF159:         END IF
169: 160: 
170:         ! Variables from common blocks or modules with intent in or inout are copied into C types161:         ! Variables from common blocks or modules with intent in or inout are copied into C types
171: 162: 
172:         DO K = 1,NRIGIDBODY163:         DO K = 1,NRIGIDBODY
173:             DO J = 1,3164:             DO J = 1,3
174:                 DO I = 1,MAXSITE165:                 DO I = 1,MAXSITE
192:         END DO183:         END DO
193: 184: 
194:         DO K = 1,3185:         DO K = 1,3
195:             DO J = 1,3186:             DO J = 1,3
196:                 DO I = 1,NRIGIDBODY187:                 DO I = 1,NRIGIDBODY
197:                    C_IINVERSE((K - 1)*3*NRIGIDBODY + (J - 1)*NRIGIDBODY + I) = IINVERSE(I,J,K)188:                    C_IINVERSE((K - 1)*3*NRIGIDBODY + (J - 1)*NRIGIDBODY + I) = IINVERSE(I,J,K)
198:                 END DO189:                 END DO
199:             END DO190:             END DO
200:         END DO191:         END DO
201: 192: 
202:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN 
203:             DO J = 1,NADDTARGET 
204:                 DO I = 1,NADDTARGET 
205:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J) 
206:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J) 
207:                 END DO 
208:             END DO 
209:         ELSE 
210:             C_LJADDREP(:) = 1.0 
211:             C_LJADDATT(:) = 1.0 
212:         END IF 
213:  
214:         DO I = 1,N193:         DO I = 1,N
215:             C_XCOORDS(I) = XCOORDS(I)194:             C_XCOORDS(I) = XCOORDS(I)
216:         END DO195:         END DO
217: 196: 
218:         DO I = 1,(N/3)197:         DO I = 1,(N/3)
219:             C_FROZEN(I) = FROZEN(I)198:             C_FROZEN(I) = FROZEN(I)
220:         END DO199:         END DO
221: 200: 
222:         C_CUDAPOT = CUDAPOT201:         C_CUDAPOT = CUDAPOT
223:         C_DEBUG = DEBUG202:         C_DEBUG = DEBUG
229:         C_NRIGIDBODY = NRIGIDBODY208:         C_NRIGIDBODY = NRIGIDBODY
230:         C_MAXSITE = MAXSITE209:         C_MAXSITE = MAXSITE
231:         C_COLDFUSIONLIMIT = COLDFUSIONLIMIT210:         C_COLDFUSIONLIMIT = COLDFUSIONLIMIT
232:         C_COLDFUSION = .FALSE.211:         C_COLDFUSION = .FALSE.
233:         C_DGUESS = DGUESS212:         C_DGUESS = DGUESS
234:         C_BQMAX = BQMAX213:         C_BQMAX = BQMAX
235:         C_FREEZE = FREEZE214:         C_FREEZE = FREEZE
236:         C_NFREEZE = NFREEZE215:         C_NFREEZE = NFREEZE
237:         PROJECT = .FALSE.216:         PROJECT = .FALSE.
238:         C_AACONVERGENCET = AACONVERGENCET217:         C_AACONVERGENCET = AACONVERGENCET
239:         C_NADDTARGET = NADDTARGET 
240: 218: 
241:         ! 'C_' prefix denotes those variables which have intent out or inout or are copies of those from common blocks/modules 219:         ! 'C_' prefix denotes those variables which have intent out or inout or are copies of those from common blocks/modules 
242:         CALL CUDA_LBFGS(N, C_XCOORDS, EPS, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, & 220:         CALL CUDA_LBFGS(N, C_XCOORDS, EPS, C_MFLAG, C_ENERGY, ITMAX, C_ITDONE, C_MAXBFGS, C_MAXERISE, C_RMS, C_CUDAPOT, & 
243:                        C_DEBUG, C_CUDATIMET, ECALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, C_DGUESS, MUPDATE, & 221:                        C_DEBUG, C_CUDATIMET, ECALLS, C_COLDFUSION, C_COLDFUSIONLIMIT, C_DGUESS, MUPDATE, & 
244:                        C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, & 222:                        C_ATOMRIGIDCOORDT, C_DEGFREEDOMS, C_NRIGIDBODY, C_NSITEPERBODY, C_RIGIDGROUPS, C_MAXSITE, & 
245:                        C_SITESRIGIDBODY, C_RIGIDSINGLES, C_BQMAX, C_IINVERSE, PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, & 223:                        C_SITESRIGIDBODY, C_RIGIDSINGLES, C_BQMAX, C_IINVERSE, PROJECT, C_FREEZE, C_FROZEN, C_NFREEZE, & 
246:                        POTENTIALTIME, C_AACONVERGENCET, C_NADDTARGET, C_LJADDREP, C_LJADDATT)224:                        POTENTIALTIME, C_AACONVERGENCET)
247: 225: 
248:         ! Make sure C types with intent out or inout are coverted back to Fortran ones226:         ! Make sure C types with intent out or inout are coverted back to Fortran ones
249: 227: 
250:         DO I = 1,N228:         DO I = 1,N
251:             XCOORDS(I) = DBLE(C_XCOORDS(I))229:             XCOORDS(I) = DBLE(C_XCOORDS(I))
252:         END DO230:         END DO
253: 231: 
254:         ENERGY = DBLE(C_ENERGY)232:         ENERGY = DBLE(C_ENERGY)
255:         RMS = DBLE(C_RMS)233:         RMS = DBLE(C_RMS)
256:         MFLAG = LOGICAL(C_MFLAG)234:         MFLAG = LOGICAL(C_MFLAG)
259: 237: 
260:         IF (COLDFUSION) THEN238:         IF (COLDFUSION) THEN
261:             WRITE(MYUNIT,'(A,G20.10)') 'ENERGY=',ENERGY239:             WRITE(MYUNIT,'(A,G20.10)') 'ENERGY=',ENERGY
262:             WRITE(MYUNIT,'(A,2G20.10)') ' Cold fusion diagnosed - step discarded; energy and threshold=',ENERGY,COLDFUSIONLIMIT240:             WRITE(MYUNIT,'(A,2G20.10)') ' Cold fusion diagnosed - step discarded; energy and threshold=',ENERGY,COLDFUSIONLIMIT
263:             ENERGY = 1.0D6241:             ENERGY = 1.0D6
264:             POTEL = 1.0D6242:             POTEL = 1.0D6
265:             RMS = 1.0D0243:             RMS = 1.0D0
266:         END IF244:         END IF
267: 245: 
268:         NPCALL = NPCALL + INT(ECALLS)246:         NPCALL = NPCALL + INT(ECALLS)
269:  
270:         IF (CENT) THEN  
271:             CALL CENTRE2(XCOORDS) 
272:         END IF 
273:  
274: #endif /* DUMMY_CUDA */247: #endif /* DUMMY_CUDA */
275: 248: 
276:     END SUBROUTINE CUDA_LBFGS_WRAPPER249:     END SUBROUTINE CUDA_LBFGS_WRAPPER
277: 250: 
278:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)251:     SUBROUTINE CUDA_ENEGRAD_WRAPPER(NATOMS, COORDS, TOTENERGY, GRADIENTS)
279: 252: 
280:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET253:         INTEGER(KIND=C_INT) :: NATOMS
281: 254: 
282:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY255:         REAL(KIND=C_DOUBLE) :: C_TOTENERGY
283:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS256:         REAL(KIND=C_DOUBLE), DIMENSION(3*NATOMS) :: COORDS, C_GRADIENTS
284: #ifndef DUMMY_CUDA257: 
285:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT258:         INTEGER :: X
286: #else ifdef DUMMY_CUDA 
287:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT 
288: #endif 
289:         INTEGER :: X, I, J 
290: 259: 
291:         DOUBLE PRECISION :: TOTENERGY260:         DOUBLE PRECISION :: TOTENERGY
292:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS261:         DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRADIENTS
293: 262: 
294: #ifndef DUMMY_CUDA263: #ifndef DUMMY_CUDA
295:  
296:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN 
297:             DO J = 1,NADDTARGET 
298:                 DO I = 1,NADDTARGET 
299:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J) 
300:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J) 
301:                 END DO 
302:             END DO 
303:         ELSE 
304:             C_LJADDREP(:) = 1.0 
305:             C_LJADDATT(:) = 1.0 
306:         END IF 
307:  
308:         C_NADDTARGET = NADDTARGET 
309:  
310:         ! Calculates the energy and gradients on the GPU using the GB potential264:         ! Calculates the energy and gradients on the GPU using the GB potential
311:         CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)265:         CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_TOTENERGY, C_GRADIENTS)
312: 266: 
313:         TOTENERGY = DBLE(C_TOTENERGY)267:         TOTENERGY = DBLE(C_TOTENERGY)
314:         268:         
315:         DO X = 1,(3*NATOMS)269:         DO X = 1,(3*NATOMS)
316:             GRADIENTS(X) = DBLE(C_GRADIENTS(X))270:             GRADIENTS(X) = DBLE(C_GRADIENTS(X))
317:         END DO271:         END DO
318: #endif272: #endif
319:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER273:     END SUBROUTINE CUDA_ENEGRAD_WRAPPER
320: 274: 
321:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)275:     SUBROUTINE CUDA_NUMERICAL_HESS(NATOMS, COORDS, HESSIAN, DELTA)
322:         IMPLICIT NONE276:         IMPLICIT NONE
323: 277: 
324:         INTEGER(KIND=C_INT) :: NATOMS, C_NADDTARGET278:         INTEGER(KIND=C_INT) :: NATOMS
325:         REAL(KIND=C_DOUBLE) :: C_ENERGY279:         REAL(KIND=C_DOUBLE) :: C_ENERGY
326:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)280:         REAL(KIND=C_DOUBLE) :: COORDS(3*NATOMS), C_GRADIENTS(3*NATOMS)
327: #ifndef DUMMY_CUDA 
328:         REAL(KIND=C_DOUBLE), DIMENSION(NADDTARGET*NADDTARGET) :: C_LJADDREP, C_LJADDATT 
329: #else ifdef DUMMY_CUDA 
330:         REAL(KIND=C_DOUBLE), DIMENSION(1, 1) :: C_LJADDREP, C_LJADDATT 
331: #endif 
332:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)281:         DOUBLE PRECISION    :: HESSIAN(3*NATOMS, 3*NATOMS)
333:         DOUBLE PRECISION    :: DELTA282:         DOUBLE PRECISION    :: DELTA
334:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)283:         DOUBLE PRECISION    :: GRAD_PLUS(3*NATOMS), GRAD_MINUS(3*NATOMS)
335:         INTEGER             :: I, J284:         INTEGER             :: I
336: 285: 
337: #ifndef DUMMY_CUDA286: #ifndef DUMMY_CUDA
338:         IF (ALLOCATED(LJADDREP) .AND. ALLOCATED(LJADDATT) .AND. LJADD3T) THEN 
339:             DO J = 1,NADDTARGET 
340:                 DO I = 1,NADDTARGET 
341:                     C_LJADDREP((J - 1)*NADDTARGET + I) = LJADDREP(I,J) 
342:                     C_LJADDATT((J - 1)*NADDTARGET + I) = LJADDATT(I,J) 
343:                 END DO 
344:             END DO 
345:         ELSE 
346:             C_LJADDREP(:) = 1.0 
347:             C_LJADDATT(:) = 1.0 
348:         END IF 
349:  
350:         C_NADDTARGET = NADDTARGET 
351:  
352:         DO I = 1, 3*NATOMS287:         DO I = 1, 3*NATOMS
353:             ! Plus288:             ! Plus
354:             COORDS(I) = COORDS(I) + DELTA289:             COORDS(I) = COORDS(I) + DELTA
355:             CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_ENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)290:             CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_ENERGY, C_GRADIENTS)
356:             GRAD_PLUS(:) = DBLE(C_GRADIENTS(:))291:             GRAD_PLUS(:) = DBLE(C_GRADIENTS(:))
357:             ! Minus292:             ! Minus
358:             COORDS(I) = COORDS(I) - 2.0D0 * DELTA293:             COORDS(I) = COORDS(I) - 2.0D0 * DELTA
359:             CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_ENERGY, C_GRADIENTS, C_NADDTARGET, C_LJADDREP, C_LJADDATT)294:             CALL CUDA_ENEGRAD_CPUTOGPU(NATOMS, COORDS, C_ENERGY, C_GRADIENTS)
360:             GRAD_MINUS(:) = DBLE(C_GRADIENTS(:))295:             GRAD_MINUS(:) = DBLE(C_GRADIENTS(:))
361:             ! Reset coords296:             ! Reset coords
362:             COORDS(I) = COORDS(I) + DELTA297:             COORDS(I) = COORDS(I) + DELTA
363:             ! Calculate hessian298:             ! Calculate hessian
364:             HESSIAN(I, :) = (GRAD_PLUS(:) - GRAD_MINUS(:)) / (2.0D0 * DELTA)299:             HESSIAN(I, :) = (GRAD_PLUS(:) - GRAD_MINUS(:)) / (2.0D0 * DELTA)
365:         END DO300:         END DO
366: #endif301: #endif
367:     END SUBROUTINE CUDA_NUMERICAL_HESS302:     END SUBROUTINE CUDA_NUMERICAL_HESS
368: 303: 
369: END MODULE MODCUDALBFGS304: END MODULE MODCUDALBFGS


r32161/mylbfgs.f90 2017-03-22 12:30:29.968265606 +0000 r32160/mylbfgs.f90 2017-03-22 12:30:33.512312362 +0000
 56:       END IF 56:       END IF
 57: ! Call to CUDA LBFGS for rigid minimisation to tolerance epsrigid before the atomistic minimisation 57: ! Call to CUDA LBFGS for rigid minimisation to tolerance epsrigid before the atomistic minimisation
 58:       IF (CUDAT) THEN 58:       IF (CUDAT) THEN
 59:          IF (.NOT. (AMBER12T)) THEN 59:          IF (.NOT. (AMBER12T)) THEN
 60:             CALL RADCOM(XCOORDS, .TRUE.) ! Evaporated atoms moved back in 60:             CALL RADCOM(XCOORDS, .TRUE.) ! Evaporated atoms moved back in
 61:          END IF 61:          END IF
 62:          CALL CUDA_LBFGS_WRAPPER(N, M, XCOORDS, EPS_TEMP, MFLAG, ENERGY, ITMAX, ITDONE, RESET) 62:          CALL CUDA_LBFGS_WRAPPER(N, M, XCOORDS, EPS_TEMP, MFLAG, ENERGY, ITMAX, ITDONE, RESET)
 63:          IF (.NOT. (AMBER12T)) THEN 63:          IF (.NOT. (AMBER12T)) THEN
 64:             EVAP = .FALSE. 64:             EVAP = .FALSE.
 65:             CALL RADCOM(XCOORDS, .FALSE.) ! Evaporated atoms NOT moved back in 65:             CALL RADCOM(XCOORDS, .FALSE.) ! Evaporated atoms NOT moved back in
 66:             IF (EVAP) THEN 66:             MFLAG = (.NOT. EVAP)
 67:                MFLAG = .FALSE. 
 68:             END IF 
 69:          END IF 67:          END IF
 70:       ELSE 68:       ELSE
 71:          !CALL MYMYLBFGS(N, M, XCOORDS, DIAGCO, EPS_TEMP, MFLAG, ENERGY, ITMAX, ITDONE, RESET, NP) !replaced by wrapper jk669  69:          !CALL MYMYLBFGS(N, M, XCOORDS, DIAGCO, EPS_TEMP, MFLAG, ENERGY, ITMAX, ITDONE, RESET, NP) !replaced by wrapper jk669 
 72:          call min_wrapper(n,m,xcoords,diagco,eps,mflag,energy,itmax,itdone,reset,np) 70:          call min_wrapper(n,m,xcoords,diagco,eps,mflag,energy,itmax,itdone,reset,np)
 73:       END IF 71:       END IF
 74:       IF (DEBUG.AND.HYBRIDMINT) WRITE(MYUNIT, '(A)') ' HYBRIDMIN> Rigid body minimisation converged, switching to all-atom' 72:       IF (DEBUG.AND.HYBRIDMINT) WRITE(MYUNIT, '(A)') ' HYBRIDMIN> Rigid body minimisation converged, switching to all-atom'
 75: ! Convert back to atomistic coordinates.  73: ! Convert back to atomistic coordinates. 
 76:       XRIGIDCOORDS(1:DEGFREEDOMS) = XCOORDS(1:DEGFREEDOMS) 74:       XRIGIDCOORDS(1:DEGFREEDOMS) = XCOORDS(1:DEGFREEDOMS)
 77:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS) 75:       CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)
 78:       ATOMRIGIDCOORDT = .TRUE. 76:       ATOMRIGIDCOORDT = .TRUE.
 83:    IF ((.NOT. RIGIDINIT) .OR. HYBRIDMINT) THEN 81:    IF ((.NOT. RIGIDINIT) .OR. HYBRIDMINT) THEN
 84:       IF (CUDAT) THEN 82:       IF (CUDAT) THEN
 85: ! Call to CUDA LBFGS for rigid minimisation 83: ! Call to CUDA LBFGS for rigid minimisation
 86:          IF (.NOT. (AMBER12T)) THEN 84:          IF (.NOT. (AMBER12T)) THEN
 87:             CALL RADCOM(XCOORDS, .TRUE.) ! Evaporated atoms moved back in 85:             CALL RADCOM(XCOORDS, .TRUE.) ! Evaporated atoms moved back in
 88:          END IF 86:          END IF
 89:          CALL CUDA_LBFGS_WRAPPER(N, M, XCOORDS, EPS, MFLAG, ENERGY, ITMAX, ITDONE, RESET) 87:          CALL CUDA_LBFGS_WRAPPER(N, M, XCOORDS, EPS, MFLAG, ENERGY, ITMAX, ITDONE, RESET)
 90:          IF (.NOT. (AMBER12T)) THEN 88:          IF (.NOT. (AMBER12T)) THEN
 91:             EVAP = .FALSE. 89:             EVAP = .FALSE.
 92:             CALL RADCOM(XCOORDS, .FALSE.) ! Evaporated atoms NOT moved back in 90:             CALL RADCOM(XCOORDS, .FALSE.) ! Evaporated atoms NOT moved back in
 93:             IF (EVAP) THEN 91:             MFLAG = (.NOT. EVAP)
 94:                MFLAG = .FALSE. 
 95:             END IF 
 96:          END IF 92:          END IF
 97:       ELSE 93:       ELSE
 98:          !CALL MYMYLBFGS(N, M, XCOORDS, DIAGCO, EPS_TEMP, MFLAG, ENERGY, ITMAX, ITDONE, RESET, NP) 94:          !CALL MYMYLBFGS(N, M, XCOORDS, DIAGCO, EPS_TEMP, MFLAG, ENERGY, ITMAX, ITDONE, RESET, NP)
 99:          call min_wrapper(n,m,xcoords,diagco,eps,mflag,energy,itmax,itdone,reset,np) !replaced by wrapper jk669  95:          call min_wrapper(n,m,xcoords,diagco,eps,mflag,energy,itmax,itdone,reset,np) !replaced by wrapper jk669 
100:       END IF 96:       END IF
101:    END IF 97:    END IF
102:  98: 
103:    GRADPROBLEMT = .FALSE. 99:    GRADPROBLEMT = .FALSE.
104:    IF (RMS.EQ.1.0D-100) THEN100:    IF (RMS.EQ.1.0D-100) THEN
105:        IF (QCIPOTT) THEN101:        IF (QCIPOTT) THEN
152:    integer(INT32), intent(inout) :: iter148:    integer(INT32), intent(inout) :: iter
153:    logical, intent(inout) :: reset149:    logical, intent(inout) :: reset
154:    integer(INT32), intent(in) :: np150:    integer(INT32), intent(in) :: np
155: 151: 
156: if (SQNMT) then !SQNM minimizer is turned on. 152: if (SQNMT) then !SQNM minimizer is turned on. 
157:    call sqnm(numcoords,xcoords,maxrmsgrad,SQNM_HISTMAX,converget,energy,itermax,iter)153:    call sqnm(numcoords,xcoords,maxrmsgrad,SQNM_HISTMAX,converget,energy,itermax,iter)
158: else154: else
159:    call MYMYLBFGS(numcoords,mupdates,xcoords,diagco,maxrmsgrad,converget,energy,itermax,iter,reset,np)155:    call MYMYLBFGS(numcoords,mupdates,xcoords,diagco,maxrmsgrad,converget,energy,itermax,iter,reset,np)
160: end if 156: end if 
161: 157: 
162: end subroutine min_wrapper158: end subroutine min_wrapper
163: 159: 


r32161/potential.f90 2017-03-22 12:30:30.188268510 +0000 r32160/potential.f90 2017-03-22 12:30:33.736315318 +0000
130:    NPCALL=NPCALL+1130:    NPCALL=NPCALL+1
131: 10 CONTINUE131: 10 CONTINUE
132: 132: 
133:    IF (MULTIPOTT) THEN133:    IF (MULTIPOTT) THEN
134:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN134:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
135:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)135:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
136:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)136:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
137:        END IF137:        END IF
138: 138: 
139: !       CALL RAD(X, GRAD, EREAL, GRADT)139: !       CALL RAD(X, GRAD, EREAL, GRADT)
140:        CALL MULTIPOT_CALL(X, GRADATOMS, EREAL, GRADT, SECT)140:        CALL MULTIPOT_CALL(X, GRAD, EREAL, GRADT, SECT)
141:        GRAD(1:3*NATOMS)=GRADATOMS(:) 
142: 141: 
143:        IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN142:        IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
144: 143: 
145:            X(DEGFREEDOMS+1:3*NATOMS)=0.0D0144:            X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
146:            X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)145:            X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
147: 146: 
148:            IF(GRADT) THEN147:            IF(GRADT) THEN
149:               CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD)148:               CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD)
150:               GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0149:               GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
151:               GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)150:               GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
183:          RETURN182:          RETURN
184:       END IF183:       END IF
185:    ELSE IF (LJADDT) THEN184:    ELSE IF (LJADDT) THEN
186: 185: 
187:       IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN186:       IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
188:           XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)187:           XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
189:           CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)188:           CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
190:       ENDIF189:       ENDIF
191: 190: 
192:       IF (LJADD3T) THEN191:       IF (LJADD3T) THEN
193:          CALL LJADD3(NATOMS, X, GRADATOMS, EREAL, GRADT, SECT)192:          CALL LJADD3(NATOMS, X, GRAD, EREAL, GRADT, SECT)
194:       ELSEIF (LJADD2T) THEN193:       ELSEIF (LJADD2T) THEN
195:          CALL LJADD2(NATOMS, X, GRADATOMS, EREAL, GRADT, SECT)194:          CALL LJADD2(NATOMS, X, GRAD, EREAL, GRADT, SECT)
196:       ELSE195:       ELSE
197:          CALL LJADD(NATOMS, X, GRADATOMS, EREAL, GRADT, SECT)196:          CALL LJADD(NATOMS, X, GRAD, EREAL, GRADT, SECT)
198:       ENDIF197:       ENDIF
199:       GRAD(1:3*NATOMS)=GRADATOMS(:) 
200: 198: 
201:       IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN199:       IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
202: 200: 
203:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0201:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
204:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)202:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
205: 203: 
206:          IF(GRADT) THEN204:          IF(GRADT) THEN
207:             CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD)205:             CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD)
208:             GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0206:             GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
209:             GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)207:             GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
497:             HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0495:             HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0
498:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)496:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) = XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)
499:          END IF497:          END IF
500:       END IF498:       END IF
501: ! AMBER 9 Energy and gradient calls499: ! AMBER 9 Energy and gradient calls
502:    ELSE IF (AMBERT) THEN500:    ELSE IF (AMBERT) THEN
503: ! hk286 > Generalised rigid body501: ! hk286 > Generalised rigid body
504: ! hk286 > If rigid coords is used, then need to convert to atom coords first,502: ! hk286 > If rigid coords is used, then need to convert to atom coords first,
505: ! hk286 > compute energies, then convert gradients & coords back to rigid503: ! hk286 > compute energies, then convert gradients & coords back to rigid
506:       IF (ATOMRIGIDCOORDT) THEN504:       IF (ATOMRIGIDCOORDT) THEN
507:          CALL AMBERENERGIES(X, GRADATOMS, EREAL, .FALSE., .FALSE.)505:          CALL AMBERENERGIES(X, GRAD, EREAL, .FALSE., .FALSE.)
508:          IF (RESTRAINLT) THEN506:          IF (RESTRAINLT) THEN
509:             CALL RESTRAINLPOTENTIAL(X, GRADATOMS, EREAL, GRADT, SECT) ! hk286507:             CALL RESTRAINLPOTENTIAL(X, GRAD, EREAL, GRADT, SECT) ! hk286
510:          END IF508:          END IF
511:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
512:       ELSE509:       ELSE
513:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)510:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
514:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)511:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)
515:          CALL AMBERENERGIES(XCOORDS, GRADATOMS, EREAL, .FALSE., .FALSE.)512:          CALL AMBERENERGIES(XCOORDS, GRADATOMS, EREAL, .FALSE., .FALSE.)
516:          IF (RESTRAINLT) THEN513:          IF (RESTRAINLT) THEN
517:             CALL RESTRAINLPOTENTIAL(XCOORDS, GRADATOMS, EREAL, GRADT, SECT) ! hk286514:             CALL RESTRAINLPOTENTIAL(XCOORDS, GRADATOMS, EREAL, GRADT, SECT) ! hk286
518:          END IF515:          END IF
519:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
520:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)516:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)
521:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)517:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
522:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0518:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
523:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)519:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
524:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0520:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
525: ! csw34> apply compression if using COMPRESSRIGID and at least one COM-COM distance521: ! csw34> apply compression if using COMPRESSRIGID and at least one COM-COM distance
526: !        is greater than RIGIDCOMDIST522: !        is greater than RIGIDCOMDIST
527:          IF (COMPRESSRIGIDT) THEN523:          IF (COMPRESSRIGIDT) THEN
528:             COMTEST=.TRUE.524:             COMTEST=.TRUE.
529:             CALL GENRIGID_DISTANCECHECK(XRIGIDCOORDS, RIGIDCOMDIST, COMTEST)525:             CALL GENRIGID_DISTANCECHECK(XRIGIDCOORDS, RIGIDCOMDIST, COMTEST)
530: ! If the threshold distance is exceeded, apply the compression526: ! If the threshold distance is exceeded, apply the compression
531:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL before compression=', EREAL527:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL before compression=', EREAL
532:             IF (.NOT.COMTEST) CALL GENRIGID_COMPRESS(XRIGIDCOORDS, GRAD, EREAL, KCOMP_RIGID)528:             IF (.NOT.COMTEST) CALL GENRIGID_COMPRESS(XRIGIDCOORDS, GRAD, EREAL, KCOMP_RIGID)
533:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL after compression=', EREAL529:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL after compression=', EREAL
 530: ! Equate GRAD and XRIGIDGRAD to allow new RMS convergence condition to work
 531:             XRIGIDGRAD(1:DEGFREEDOMS)=GRAD(1:DEGFREEDOMS)
534:          END IF532:          END IF
535:       END IF533:       END IF
536: 534: 
537:       IF (SECT) THEN535:       IF (SECT) THEN
538:          VNEW=0.0D0536:          VNEW=0.0D0
539: ! khs26> Copied analytical second derivatives from OPTIM537: ! khs26> Copied analytical second derivatives from OPTIM
540:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN538:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
541:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)539:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
542:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)540:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
543:          END IF541:          END IF
607:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) =XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)605:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) =XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS)
608:          END IF606:          END IF
609:       END IF607:       END IF
610: 608: 
611: ! hk286609: ! hk286
612:    ELSE IF (CHRMMT) THEN610:    ELSE IF (CHRMMT) THEN
613: ! hk286 > Generalised rigid body611: ! hk286 > Generalised rigid body
614: ! hk286 > If rigid coords is used, then need to convert to atom coords first,612: ! hk286 > If rigid coords is used, then need to convert to atom coords first,
615: ! hk286 > compute energies, then convert gradients & coords back to rigid613: ! hk286 > compute energies, then convert gradients & coords back to rigid
616:       IF (ATOMRIGIDCOORDT) THEN614:       IF (ATOMRIGIDCOORDT) THEN
617:          CALL OCHARMM(X, GRADATOMS, EREAL, GRADT)615:          CALL OCHARMM(X, GRAD, EREAL, GRADT)
618:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
619:       ELSE616:       ELSE
620:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)617:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
621:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)618:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)
622:          CALL OCHARMM(XCOORDS, GRADATOMS, EREAL, GRADT)619:          CALL OCHARMM(XCOORDS, GRADATOMS, EREAL, GRADT)
623:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
624:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)620:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)
625:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)621:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
626:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0622:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
627:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)623:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
628:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0624:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
629: ! csw34> apply compression if using COMPRESSRIGID and at least one COM-COM distance625: ! csw34> apply compression if using COMPRESSRIGID and at least one COM-COM distance
630: !        is greater than RIGIDCOMDIST626: !        is greater than RIGIDCOMDIST
631:          IF (COMPRESSRIGIDT) THEN627:          IF (COMPRESSRIGIDT) THEN
632:             COMTEST=.TRUE.628:             COMTEST=.TRUE.
633:             CALL GENRIGID_DISTANCECHECK(XRIGIDCOORDS, RIGIDCOMDIST, COMTEST)629:             CALL GENRIGID_DISTANCECHECK(XRIGIDCOORDS, RIGIDCOMDIST, COMTEST)
634: ! If the threshold distance is exceeded, apply the compression630: ! If the threshold distance is exceeded, apply the compression
635:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL before compression=', EREAL631:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL before compression=', EREAL
636:             IF (.NOT.COMTEST) CALL GENRIGID_COMPRESS(XRIGIDCOORDS, GRAD, EREAL, KCOMP_RIGID)632:             IF (.NOT.COMTEST) CALL GENRIGID_COMPRESS(XRIGIDCOORDS, GRAD, EREAL, KCOMP_RIGID)
637:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL after compression=', EREAL633:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL after compression=', EREAL
 634: ! Equate GRAD and XRIGIDGRAD to allow new RMS convergence condition to work
 635:             XRIGIDGRAD(1:DEGFREEDOMS)=GRAD(1:DEGFREEDOMS)
638:          END IF636:          END IF
639:       END IF637:       END IF
640: 638: 
641: ! hk286639: ! hk286
642:    ELSE IF (DMACRYST) THEN640:    ELSE IF (DMACRYST) THEN
643:       IF (.NOT. ATOMRIGIDCOORDT) THEN641:       IF (.NOT. ATOMRIGIDCOORDT) THEN
644:          CALL DMACRYS_POTENTIAL(X, GRAD, EREAL, GRADT)642:          CALL DMACRYS_POTENTIAL(X, GRAD, EREAL, GRADT)
645:       ELSE643:       ELSE
646:          CALL TRANSFORMCTORIGID(X, XRIGIDCOORDS)644:          CALL TRANSFORMCTORIGID(X, XRIGIDCOORDS)
647:          CALL DMACRYS_POTENTIAL(XRIGIDCOORDS, GRAD, EREAL, GRADT)645:          CALL DMACRYS_POTENTIAL(XRIGIDCOORDS, GRAD, EREAL, GRADT)
650:    ELSE IF (DZTEST) THEN648:    ELSE IF (DZTEST) THEN
651:       CALL RAD(X, GRAD, EREAL, GRADT)649:       CALL RAD(X, GRAD, EREAL, GRADT)
652:       CALL DZPOT(X, GRAD, EREAL, GRADT, SECT)650:       CALL DZPOT(X, GRAD, EREAL, GRADT, SECT)
653: ! vr274: userpot, linkage dependant potential651: ! vr274: userpot, linkage dependant potential
654: 652: 
655:    ELSE IF (USERPOTT) THEN653:    ELSE IF (USERPOTT) THEN
656: ! hk286 > Generalised rigid body654: ! hk286 > Generalised rigid body
657: ! hk286 > If rigid coords is used, then need to convert to atom coords first,655: ! hk286 > If rigid coords is used, then need to convert to atom coords first,
658: ! hk286 > compute energies, then convert gradients & coords back to rigid656: ! hk286 > compute energies, then convert gradients & coords back to rigid
659:       IF (ATOMRIGIDCOORDT) THEN657:       IF (ATOMRIGIDCOORDT) THEN
660:          CALL USERPOT_POTENTIAL(3*NATOMS, X, GRADATOMS, EREAL, GRADT)658:          CALL USERPOT_POTENTIAL(3*NATOMS, X, GRAD, EREAL, GRADT)
661:          IF (RESTRAINLT) THEN659:          IF (RESTRAINLT) THEN
662:             CALL RESTRAINLPOTENTIAL(X, GRADATOMS, EREAL, GRADT, SECT) ! hk286660:             CALL RESTRAINLPOTENTIAL(X, GRAD, EREAL, GRADT, SECT) ! hk286
663:          END IF661:          END IF
664:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
665:       ELSE662:       ELSE
666:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)663:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
667:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)664:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)
668:          CALL USERPOT_POTENTIAL(3*NATOMS, XCOORDS, GRADATOMS, EREAL, GRADT)665:          CALL USERPOT_POTENTIAL(3*NATOMS, XCOORDS, GRADATOMS, EREAL, GRADT)
669:          IF (RESTRAINLT) THEN666:          IF (RESTRAINLT) THEN
670:             CALL RESTRAINLPOTENTIAL(XCOORDS, GRADATOMS, EREAL, GRADT, SECT) ! hk286667:             CALL RESTRAINLPOTENTIAL(XCOORDS, GRADATOMS, EREAL, GRADT, SECT) ! hk286
671:          END IF668:          END IF
672:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
673:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)669:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)
674:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)670:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
675:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0671:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
676:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)672:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
677:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0673:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
678: ! csw34> apply compression if using COMPRESSRIGID and at least one COM-COM distance674: ! csw34> apply compression if using COMPRESSRIGID and at least one COM-COM distance
679: !        is greater than RIGIDCOMDIST675: !        is greater than RIGIDCOMDIST
680:          IF (COMPRESSRIGIDT) THEN676:          IF (COMPRESSRIGIDT) THEN
681:             COMTEST=.TRUE.677:             COMTEST=.TRUE.
682:             CALL GENRIGID_DISTANCECHECK(XRIGIDCOORDS, RIGIDCOMDIST, COMTEST)678:             CALL GENRIGID_DISTANCECHECK(XRIGIDCOORDS, RIGIDCOMDIST, COMTEST)
683: ! If the threshold distance is exceeded, apply the compression679: ! If the threshold distance is exceeded, apply the compression
684:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL before compression=', EREAL680:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL before compression=', EREAL
685:             IF (.NOT.COMTEST) CALL GENRIGID_COMPRESS(XRIGIDCOORDS, GRAD, EREAL, KCOMP_RIGID)681:             IF (.NOT.COMTEST) CALL GENRIGID_COMPRESS(XRIGIDCOORDS, GRAD, EREAL, KCOMP_RIGID)
686:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL after compression=', EREAL682:             IF (DEBUG) WRITE(MYUNIT,'(A, F20.10)' ) 'EREAL after compression=', EREAL
 683: ! Equate GRAD and XRIGIDGRAD to allow new RMS convergence condition to work
 684:             XRIGIDGRAD(1:DEGFREEDOMS)=GRAD(1:DEGFREEDOMS)
687:          END IF685:          END IF
688:       END IF686:       END IF
689: 687: 
690:       IF (SECT) THEN688:       IF (SECT) THEN
691:          VNEW=0.0D0689:          VNEW=0.0D0
692: ! khs26> Copied analytical second derivatives from OPTIM690: ! khs26> Copied analytical second derivatives from OPTIM
693:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN691:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
694:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)692:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
695:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)693:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
696:          END IF694:          END IF
756:    ELSE IF (BGUPTAT) THEN754:    ELSE IF (BGUPTAT) THEN
757:       CALL RAD(X, GRAD, EREAL, GRADT)755:       CALL RAD(X, GRAD, EREAL, GRADT)
758:       CALL BGUPTA(X, GRAD, EREAL, GRADT)756:       CALL BGUPTA(X, GRAD, EREAL, GRADT)
759: 757: 
760:    ELSE IF (DJWRBT) THEN758:    ELSE IF (DJWRBT) THEN
761:       IF (DJWRBID.EQ.1) THEN759:       IF (DJWRBID.EQ.1) THEN
762:          IF (.NOT. ATOMRIGIDCOORDT) THEN760:          IF (.NOT. ATOMRIGIDCOORDT) THEN
763:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)761:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
764:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)762:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)
765:          ENDIF763:          ENDIF
766:          CALL DJWGR1(NATOMS,X,GRADATOMS,EREAL,GRADT)764:          CALL DJWGR1(NATOMS,X,GRAD,EREAL,GRADT)
767:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
768:          IF (.NOT.ATOMRIGIDCOORDT) THEN765:          IF (.NOT.ATOMRIGIDCOORDT) THEN
769:              X(DEGFREEDOMS+1:3*NATOMS)=0.0D0766:              X(DEGFREEDOMS+1:3*NATOMS)=0.0D0
770:              X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)767:              X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS)
771:              IF (GRADT) THEN768:              IF (GRADT) THEN
772:                 CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD)769:                 CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD)
773:                 GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0770:                 GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0
774:                 GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)771:                 GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
775:              ENDIF772:              ENDIF
776:          ENDIF773:          ENDIF
777:       ENDIF774:       ENDIF
1424: !            J2=3*J11421: !            J2=3*J1
1425: !            GRAD(J2-2)=GRAD(J2-2)-XG1422: !            GRAD(J2-2)=GRAD(J2-2)-XG
1426: !            GRAD(J2-1)=GRAD(J2-1)-YG1423: !            GRAD(J2-1)=GRAD(J2-1)-YG
1427: !            GRAD(J2)=  GRAD(J2)-ZG1424: !            GRAD(J2)=  GRAD(J2)-ZG
1428: !         END DO1425: !         END DO
1429: !      END IF1426: !      END IF
1430: 1427: 
1431:       IF (CSMT .AND. (.NOT.SYMMETRIZECSM)) THEN1428:       IF (CSMT .AND. (.NOT.SYMMETRIZECSM)) THEN
1432:          DUMMY2=0.0D01429:          DUMMY2=0.0D0
1433:          RMS=0.0D01430:          RMS=0.0D0
1434:       ELSE IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN1431:       ELSE IF (AACONVERGENCET .AND. (ATOMRIGIDCOORDT)) THEN
1435:          DUMMY2=SUM(GRAD(1:DEGFREEDOMS)**2)1432:          DUMMY2=SUM(GRAD(1:DEGFREEDOMS)**2)
1436:          RMS=MAX(SQRT(DUMMY2/DEGFREEDOMS), 1.0D-100)1433:          RMS=MAX(SQRT(DUMMY2/DEGFREEDOMS), 1.0D-100)
1437:          IF ((RMS < 5.0D0 * BQMAX) .AND. AACONVERGENCET .AND. (.NOT. COMPRESSRIGIDT) .AND. (.NOT. COMPON)) THEN1434:          IF (RMS < 5.0D0 * BQMAX) THEN
1438:             CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, RMS)1435:             CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, RMS)
1439:          END IF1436:          END IF
1440:       ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT) THEN1437:       ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT) THEN
1441:          DUMMY2=SUM(GRAD(1:NATOMS)**2)1438:          DUMMY2=SUM(GRAD(1:NATOMS)**2)
1442:          RMS=MAX(DSQRT(DUMMY2/(NATOMS)), 1.0D-100)1439:          RMS=MAX(DSQRT(DUMMY2/(NATOMS)), 1.0D-100)
1443:       ELSE IF (.NOT.THOMSONT) THEN1440:       ELSE IF (.NOT.THOMSONT) THEN
1444:          DUMMY2=SUM(GRAD(1:3*NATOMS)**2)1441:          DUMMY2=SUM(GRAD(1:3*NATOMS)**2)
1445:          RMS=MAX(DSQRT(DUMMY2/(3*NATOMS)), 1.0D-100)1442:          RMS=MAX(DSQRT(DUMMY2/(3*NATOMS)), 1.0D-100)
1446:       ELSE1443:       ELSE
1447:          DUMMY2=SUM(GRAD(1:2*NATOMS)**2)1444:          DUMMY2=SUM(GRAD(1:2*NATOMS)**2)


r32161/potential.h 2017-03-22 12:30:28.644248139 +0000 r32160/potential.h 2017-03-22 12:30:32.180294789 +0000
 10: #ifndef POTENTIAL_H 10: #ifndef POTENTIAL_H
 11: #define POTENTIAL_H 11: #define POTENTIAL_H
 12:  12: 
 13: class LjPotential : public CostFunction 13: class LjPotential : public CostFunction
 14: { 14: {
 15:         public: 15:         public:
 16:                 LjPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom,  16:                 LjPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom, 
 17:                                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody,  17:                                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody, 
 18:                                 int *rigidSingles, double *rigidInverse, double *coords, int nSecDiag, bool isAtomisticNotRigid,  18:                                 int *rigidSingles, double *rigidInverse, double *coords, int nSecDiag, bool isAtomisticNotRigid, 
 19:                                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze,  19:                                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze, 
 20:                                 bool isAaConvergence, int nAddTarget, double *ljAddRep, double *ljAddAtt); 20:                                 bool isAaConvergence);
 21:  21: 
 22:                 void computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf); 22:                 void computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf);
 23: }; 23: };
 24:  24: 
 25: class AmberPotential : public CostFunction 25: class AmberPotential : public CostFunction
 26: { 26: {
 27:         public: 27:         public:
 28:                 AmberPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom,  28:                 AmberPotential(Printing &debugPrinting, Timer &timer_potential, Cublas &cublas, size_t numDimensions, int nDegFreedom, 
 29:                                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody,  29:                                 int nRigidBody, int rigidMaxSite, int *nRigidSitesPerBody, int *rigidGroups, double *sitesRigidBody, 
 30:                                 int *rigidSingles, double *rigidInverse, double *x, int nSecDiag, bool isAtomisticNotRigid,  30:                                 int *rigidSingles, double *rigidInverse, double *x, int nSecDiag, bool isAtomisticNotRigid, 
 31:                                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze,  31:                                 double aaConvThreshold, double coldFusionLim, bool shouldFreeze, bool *isAtomFrozen, int nFreeze, 
 32:                                 bool isAaConvergence, int nAddTarget, double *ljAddRep, double *ljAddAtt); 32:                                 bool isAaConvergence);
 33:  33: 
 34:                 void computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf); 34:                 void computeEnergyAndGradient(const double *d_x, double *d_f, double *d_gradf);
 35: }; 35: };
 36:  36: 
 37: #endif /* end of include guard: POTENTIAL_H */ 37: #endif /* end of include guard: POTENTIAL_H */


r32161/rigid_bodies.cu 2017-03-22 12:30:28.204242333 +0000 r32160/rigid_bodies.cu 2017-03-22 12:30:31.736288932 +0000
  3:  * File rigid_bodies.cu: Implementation of rigid body transformation functions from class CostFunction.  3:  * File rigid_bodies.cu: Implementation of rigid body transformation functions from class CostFunction.
  4:  *  4:  *
  5:  **/  5:  **/
  6:   6: 
  7: #include "cost_function.h"  7: #include "cost_function.h"
  8:   8: 
  9: namespace gpu_rigid_bodies  9: namespace gpu_rigid_bodies
 10: { 10: {
 11:         // Variables on the GPU.  11:         // Variables on the GPU. 
 12:  12: 
  13:         __device__ int  singlesThreads;
  14:         __device__ int  thisBody;
  15:         __device__ int  arraySize;
  16:         __device__ int  gridSize;
  17:         __device__ int  roundMaxSite;
  18:         __device__ int  outSize;
  19:         __device__ int  nBlocks;
 13:         __device__ bool shouldFindDeriv; 20:         __device__ bool shouldFindDeriv;
 14:  21: 
 15:  22: 
 16:         // Kernels.  23:         // Kernels. 
 17:  24: 
 18:         __global__ void transform(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi); 25:         __global__ void transform(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi);
 19:  26: 
 20:         __global__ void transform2(double *d_x, const double *m_d_xRigid, const int *m_d_nRigidBody,  27:         __global__ void transform2(double *d_x, const double *m_d_xRigid, const int *m_d_nRigidBody, 
 21:                         const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  28:                         const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
 22:                         const double *m_d_sitesRigidBody, const double *d_grmi, const int *m_d_rigidMaxSite); 29:                         const double *m_d_sitesRigidBody, const double *d_grmi, const int *m_d_rigidMaxSite);
 23:  30: 
 24:         __global__ void singleAtoms(double *d_x, const double *m_d_xRigid, const int *m_d_rigidSingles, const int *m_d_nRigidBody,  31:         __global__ void singleAtoms(double *d_x, const double *m_d_xRigid, const int *m_d_rigidSingles, const int *m_d_nRigidBody);
 25:                         const int singlesThreads); 
 26:  32: 
 27:         __global__ void gradTransform1a(const double *m_d_xRigid, double *d_grmi1, double *d_grmi2, double *d_grmi3,  33:         __global__ void gradTransform1a(const double *m_d_xRigid, double *d_grmi1, double *d_grmi2, double *d_grmi3, 
 28:                         const int *m_d_nRigidBody); 34:                         const int *m_d_nRigidBody);
 29:  35: 
 30:         __global__ void gradTransform1b(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi0); 36:         __global__ void gradTransform1b(const int *m_d_nRigidBody, const double *m_d_xRigid, double *d_grmi0);
 31:  37: 
 32:         __global__ void intermediate(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi1,  38:         __global__ void intermediate(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi1, 
 33:                         const double *d_grmi2, const double *d_grmi3, const double *m_d_sitesRigidBody,  39:                         const double *d_grmi2, const double *d_grmi3, const double *m_d_sitesRigidBody, 
 34:                         const int *m_d_rigidGroups, double *d_tempArray, const int *m_d_nRigidBody,  40:                         const int *m_d_rigidGroups, double *d_tempArray, const int *m_d_nRigidBody, 
 35:                         const int *m_d_rigidMaxSite); 41:                         const int *m_d_rigidMaxSite);
 36:  42: 
 37:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles,  43:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles, 
 38:                         const int *m_d_nRigidBody, const int singlesThreads); 44:                         const int *m_d_nRigidBody);
 39:  45: 
 40:         __global__ void aaConvRmdrvt(double *d_rmi10, double *d_rmi20, double *d_rmi30); 46:         __global__ void aaConvRmdrvt(double *d_rmi10, double *d_rmi20, double *d_rmi30);
 41:  47: 
 42:         __global__ void intermediate2(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi0,  48:         __global__ void intermediate2(const double *d_gk, const int *m_d_nRigidSitesPerBody, const double *d_grmi0, 
 43:                         const double *d_grmi10, const double *d_grmi20, const double *d_grmi30,  49:                         const double *d_grmi10, const double *d_grmi20, const double *d_grmi30, 
 44:                         const double *m_d_sitesRigidBody, const int *m_d_rigidGroups, double *d_tempArray,  50:                         const double *m_d_sitesRigidBody, const int *m_d_rigidGroups, double *d_tempArray, 
 45:                         const int *m_d_nRigidBody, const int *m_d_rigidMaxSite); 51:                         const int *m_d_nRigidBody, const int *m_d_rigidMaxSite);
 46:  52: 
 47:         __global__ void aaConvTorque(const double *d_torques, const double *m_d_gkRigid, const double *d_grmi0,  53:         __global__ void aaConvTorque(const double *d_torques, const double *m_d_gkRigid, const double *d_grmi0, 
 48:                         const int *m_d_nRigidSitesPerBody, const double *m_d_rigidInverse, double *d_rmsArray,  54:                         const int *m_d_nRigidSitesPerBody, const double *m_d_rigidInverse, double *d_rmsArray, 
 51:         __global__ void warpReduce1(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  57:         __global__ void warpReduce1(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
 52:                         const int *m_d_rigidMaxSite, const double *d_gk, double *m_d_gkRigid); 58:                         const int *m_d_rigidMaxSite, const double *d_gk, double *m_d_gkRigid);
 53:  59: 
 54:         __global__ void warpReduce2(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite,  60:         __global__ void warpReduce2(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite, 
 55:                         double *m_d_gkRigid, const double *d_tempArray); 61:                         double *m_d_gkRigid, const double *d_tempArray);
 56:  62: 
 57:         __global__ void warpReduce3(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite,  63:         __global__ void warpReduce3(const int *m_d_nRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidMaxSite, 
 58:                         double *d_torques, const double *d_tempArray); 64:                         double *d_torques, const double *d_tempArray);
 59:  65: 
 60:         __global__ void fullReduce1a(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,  66:         __global__ void fullReduce1a(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
 61:                         const int *m_d_rigidMaxSite, const double *d_gk, const int *m_d_largeRigidIndices, double *d_outArray,  67:                         const int *m_d_rigidMaxSite, const double *d_gk, const int *m_d_largeRigidIndices, double3 *d_outArray);
 62:                         const int roundedMaxSite); 
 63:  68: 
 64:         __global__ void fullReduce2a(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices,  69:         __global__ void fullReduce2a(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, 
 65:                         double *d_outArray, const int roundedMaxSite, const int outSize, const int blocks); 70:                         double3 *d_outArray);
 66:  71: 
 67:         __global__ void fullReduce1b(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 72:         __global__ void fullReduce1b(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups,
 68:                         const int *m_d_rigidMaxSite, const int *m_d_largeRigidIndices, double *d_outArray, 73:                         const int *m_d_rigidMaxSite, const int *m_d_largeRigidIndices, double3 *d_outArray,
 69:                         const double *d_tempArray, const int roundedMaxSite); 74:                         const double *d_tempArray);
 70:  75: 
 71:         __global__ void fullReduce2b(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices,  76:         __global__ void fullReduce2b(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, 
 72:                         double *d_outArray, const int *m_d_nRigidBody, const int roundedMaxSite, const int outSize,  77:                         double3 *d_outArray, const int *m_d_nRigidBody);
 73:                         const int blocks); 
 74:  78: 
 75:         __global__ void fullReduce2c(const int *m_d_nLargeRigidBody, double *torques, const int *m_d_largeRigidIndices,  79:         __global__ void fullReduce2c(const int *m_d_nLargeRigidBody, double *torques, const int *m_d_largeRigidIndices, 
 76:                         double *d_outArray, const int *m_d_nRigidBody, const int roundedMaxSite, const int outSize, const int blocks); 80:                         double3 *d_outArray, const int *m_d_nRigidBody);
 77:  81: 
 78:         __global__ void fullReduce(const double *in, double* out, const int N); 82:         __global__ void fullReduce(double *in, double* out, int N);
 79:  83: 
 80:         __global__ void fullDotReduce(const double *in, double* out, const int N); 84:         __global__ void fullDotReduce(double *in, double* out, int N);
 81:  85: 
 82: } 86: }
 83:  87: 
 84:  88: 
 85:  89: 
 86: void CostFunction::transformRigidToC(double *d_x) 90: void CostFunction::transformRigidToC(double *d_x)
 87: { 91: {
 88:         using namespace gpu_rigid_bodies; 92:         using namespace gpu_rigid_bodies;
 89:  93: 
 90:         // Copy atomistic coordinates to rigid coordinates.  94:         // Copy atomistic coordinates to rigid coordinates. 
115: 119: 
116:         // Rigid body to atomistic mapping. 120:         // Rigid body to atomistic mapping. 
117:         gpu_rigid_bodies::transform2<<<gridDim, blockDim>>>(d_x, m_d_xRigid, m_d_nRigidBody, m_d_nRigidSitesPerBody, 121:         gpu_rigid_bodies::transform2<<<gridDim, blockDim>>>(d_x, m_d_xRigid, m_d_nRigidBody, m_d_nRigidSitesPerBody, 
118:                         m_d_rigidGroups, m_d_sitesRigidBody, d_grmi, m_d_rigidMaxSite); 122:                         m_d_rigidGroups, m_d_sitesRigidBody, d_grmi, m_d_rigidMaxSite); 
119:         CudaCheckError();123:         CudaCheckError();
120:         cudaDeviceSynchronize();124:         cudaDeviceSynchronize();
121: 125: 
122:         CudaSafeCall( cudaFree(d_grmi) );126:         CudaSafeCall( cudaFree(d_grmi) );
123: 127: 
124:         if (m_nDegFreedom > 6*m_nRigidBody) {128:         if (m_nDegFreedom > 6*m_nRigidBody) {
125:                 int singlesThreads = (m_nDegFreedom - 6*m_nRigidBody)/3;129:                 int hostSinglesThreads = (m_nDegFreedom - 6*m_nRigidBody)/3;
 130:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::singlesThreads, &hostSinglesThreads,  sizeof(int)) );
126: 131: 
127:                 blockDim.x = 1024;132:                 blockDim.x = 1024;
128:                 gridDim.x = (singlesThreads + blockDim.x - 1)/blockDim.x;133:                 gridDim.x = (hostSinglesThreads + blockDim.x - 1)/blockDim.x;
129: 134: 
130:                 // Treatment of single atoms not in rigid bodies. 135:                 // Treatment of single atoms not in rigid bodies. 
131:                 gpu_rigid_bodies::singleAtoms<<<gridDim, blockDim>>>(d_x, m_d_xRigid, m_d_rigidSingles, m_d_nRigidBody, singlesThreads);136:                 gpu_rigid_bodies::singleAtoms<<<gridDim, blockDim>>>(d_x, m_d_xRigid, m_d_rigidSingles, m_d_nRigidBody);
132:                 CudaCheckError();137:                 CudaCheckError();
133:                 cudaDeviceSynchronize();138:                 cudaDeviceSynchronize();
134:         }139:         }
135: }140: }
136: 141: 
137: 142: 
138: 143: 
139: void CostFunction::transformGrad(double *d_gk, double *d_x)144: void CostFunction::transformGrad(double *d_gk, double *d_x)
140: {145: {
141:         using namespace gpu_rigid_bodies;146:         using namespace gpu_rigid_bodies;
182:         CudaCheckError();187:         CudaCheckError();
183:         cudaDeviceSynchronize();188:         cudaDeviceSynchronize();
184: 189: 
185:         // Reduction for rigid bodies with more than 32 sites. 190:         // Reduction for rigid bodies with more than 32 sites. 
186:         if (m_nLargeRigidBody != 0){191:         if (m_nLargeRigidBody != 0){
187:                 int blockSize = 1024;192:                 int blockSize = 1024;
188:                 blockDim.x = blockSize;193:                 blockDim.x = blockSize;
189: 194: 
190:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies. 195:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies. 
191:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize;196:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize;
 197:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) );
192:                 int numThreads = roundedMaxSite*m_nLargeRigidBody;198:                 int numThreads = roundedMaxSite*m_nLargeRigidBody;
193: 199: 
194:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x;200:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x;
195:                 gridDim.x = blocks;201:                 gridDim.x = blocks;
196: 202: 
197:                 double *d_outArray;203:                 double3 *d_outArray;
198:                 CudaSafeCall( cudaMalloc(&d_outArray, 3 * blocks * sizeof(double)) );204:                 CudaSafeCall( cudaMalloc(&d_outArray, blocks * sizeof(double)) );
199: 205: 
200:                 gpu_rigid_bodies::fullReduce1a<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, 206:                 gpu_rigid_bodies::fullReduce1a<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, 
201:                                 m_d_rigidGroups, m_d_rigidMaxSite, d_gk, m_d_largeRigidIndices, d_outArray, roundedMaxSite);207:                                 m_d_rigidGroups, m_d_rigidMaxSite, d_gk, m_d_largeRigidIndices, d_outArray);
202:                 CudaCheckError(); 
203:                 cudaDeviceSynchronize(); 
204: 208: 
205:                 while (blocks > m_nLargeRigidBody) {209:                 while (blocks > m_nLargeRigidBody) {
206:                         int outSize = blocks;210:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::outSize, &blocks, sizeof(int)) );
207: 211: 
208:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize;212:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize;
 213:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) );
209:                         numThreads = roundedMaxSite*m_nLargeRigidBody;214:                         numThreads = roundedMaxSite*m_nLargeRigidBody;
210: 215: 
211:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x;216:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x;
212:                         gridDim.x = blocks;217:                         gridDim.x = blocks;
 218:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::nBlocks, &blocks, sizeof(int)) );
213: 219: 
214:                         gpu_rigid_bodies::fullReduce2a<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_gkRigid, 220:                         gpu_rigid_bodies::fullReduce2a<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_gkRigid, 
215:                                         m_d_largeRigidIndices, d_outArray, roundedMaxSite, outSize, blocks);221:                                         m_d_largeRigidIndices, d_outArray);
216:                         CudaCheckError(); 
217:                         cudaDeviceSynchronize(); 
218:                 }222:                 }
219: 223: 
220:                 CudaSafeCall( cudaFree(d_outArray) );224:                 CudaSafeCall( cudaFree(d_outArray) );
221:         }225:         }
222: 226: 
223:         // Temporary global memory storage for atomistic components of rigid body forces. 227:         // Temporary global memory storage for atomistic components of rigid body forces. 
224:         double *d_tempArray;228:         double *d_tempArray;
225:         int tempArraySize = 3 * m_nRigidBody * m_rigidMaxSite;229:         int tempArraySize = 3 * m_nRigidBody * m_rigidMaxSite;
226:         CudaSafeCall( cudaMalloc(&d_tempArray, tempArraySize * sizeof(double)) );230:         CudaSafeCall( cudaMalloc(&d_tempArray, tempArraySize * sizeof(double)) );
227: 231: 
251:         CudaCheckError();255:         CudaCheckError();
252:         cudaDeviceSynchronize();256:         cudaDeviceSynchronize();
253: 257: 
254:         // Reduction for rigid bodies with more than 32 sites.258:         // Reduction for rigid bodies with more than 32 sites.
255:         if (m_nLargeRigidBody != 0){259:         if (m_nLargeRigidBody != 0){
256:                 int blockSize = 1024;260:                 int blockSize = 1024;
257:                 blockDim.x = blockSize;261:                 blockDim.x = blockSize;
258: 262: 
259:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies.263:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies.
260:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize;264:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize;
 265:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) );
261:                 int numThreads = roundedMaxSite*m_nLargeRigidBody;266:                 int numThreads = roundedMaxSite*m_nLargeRigidBody;
262: 267: 
263:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x;268:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x;
264:                 gridDim.x = blocks;269:                 gridDim.x = blocks;
265: 270: 
266:                 double *d_outArray;271:                 double3 *d_outArray;
267:                 CudaSafeCall( cudaMalloc(&d_outArray, 3 * blocks * sizeof(double)) );272:                 CudaSafeCall( cudaMalloc(&d_outArray, blocks * sizeof(double)) );
268: 273: 
269:                 gpu_rigid_bodies::fullReduce1b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups,274:                 gpu_rigid_bodies::fullReduce1b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups,
270:                                 m_d_rigidMaxSite, m_d_largeRigidIndices, d_outArray, d_tempArray, roundedMaxSite);275:                                 m_d_rigidMaxSite, m_d_largeRigidIndices, d_outArray, d_tempArray);
271:                 CudaCheckError(); 
272:                 cudaDeviceSynchronize(); 
273: 276: 
274:                 while (blocks > m_nLargeRigidBody) {277:                 while (blocks > m_nLargeRigidBody) {
275:                         int outSize = blocks;278:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::outSize, &blocks, sizeof(int)) );
276: 279: 
277:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize;280:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize;
 281:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) );
278:                         numThreads = roundedMaxSite*m_nLargeRigidBody;282:                         numThreads = roundedMaxSite*m_nLargeRigidBody;
279: 283: 
280:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x;284:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x;
281:                         gridDim.x = blocks;285:                         gridDim.x = blocks;
 286:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::nBlocks, &blocks, sizeof(int)) );
282: 287: 
283:                         gpu_rigid_bodies::fullReduce2b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_gkRigid, 288:                         gpu_rigid_bodies::fullReduce2b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_gkRigid, m_d_largeRigidIndices, 
284:                                         m_d_largeRigidIndices, d_outArray, m_d_nRigidBody, roundedMaxSite, outSize, blocks);289:                                         d_outArray, m_d_nRigidBody);
285:                         CudaCheckError(); 
286:                         cudaDeviceSynchronize(); 
287:                 }290:                 }
288: 291: 
289:                 CudaSafeCall( cudaFree(d_outArray) );292:                 CudaSafeCall( cudaFree(d_outArray) );
290:         }293:         }
291: 294: 
292:         CudaSafeCall( cudaFree(d_tempArray) );295:         CudaSafeCall( cudaFree(d_tempArray) );
293: 296: 
294:         if (m_nDegFreedom > 6*m_nRigidBody) {297:         if (m_nDegFreedom > 6*m_nRigidBody) {
295:                 int singlesThreads = (m_nDegFreedom - 6*m_nRigidBody)/3;298:                 int hostSinglesThreads = (m_nDegFreedom - 6*m_nRigidBody)/3;
 299:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::singlesThreads, &hostSinglesThreads, sizeof(int)) );
296: 300: 
297:                 blockDim.x = 1024;301:                 blockDim.x = 1024;
298:                 gridDim.x = (singlesThreads + blockDim.x - 1)/blockDim.x;302:                 gridDim.x = (hostSinglesThreads + blockDim.x - 1)/blockDim.x;
299: 303: 
300:                 // Treatment of single atoms not in rigid bodies. 304:                 // Treatment of single atoms not in rigid bodies. 
301:                 gpu_rigid_bodies::gradSingleAtoms<<<gridDim, blockDim>>>(d_gk, m_d_gkRigid, m_d_rigidSingles, 305:                 gpu_rigid_bodies::gradSingleAtoms<<<gridDim, blockDim>>>(d_gk, m_d_gkRigid, m_d_rigidSingles, 
302:                                 m_d_nRigidBody, singlesThreads);306:                                 m_d_nRigidBody);
303:                 CudaCheckError();307:                 CudaCheckError();
304:                 cudaDeviceSynchronize();308:                 cudaDeviceSynchronize();
305:         }309:         }
306: 310: 
307:         // Copy rigid coordinates/gradient to atomistic coordinates/gradient and pad with zeros. 311:         // Copy rigid coordinates/gradient to atomistic coordinates/gradient and pad with zeros. 
308:         CudaSafeCall( cudaMemcpy(d_x,  zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );312:         CudaSafeCall( cudaMemcpy(d_x,  zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
309:         CudaSafeCall( cudaMemcpy(d_gk, zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );313:         CudaSafeCall( cudaMemcpy(d_gk, zeros, numDimensions * sizeof(double), cudaMemcpyHostToDevice) );
310: 314: 
311:         CudaSafeCall( cudaMemcpy(d_x,  m_d_xRigid,  m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );315:         CudaSafeCall( cudaMemcpy(d_x,  m_d_xRigid,  m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );
312:         CudaSafeCall( cudaMemcpy(d_gk, m_d_gkRigid, m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );316:         CudaSafeCall( cudaMemcpy(d_gk, m_d_gkRigid, m_nDegFreedom * sizeof(double), cudaMemcpyDeviceToDevice) );
313: 317: 
314:         delete [] zeros;318:         delete [] zeros;
 319: 
315: }320: }
316: 321: 
317: 322: 
318: 323: 
319: void CostFunction::aaConvergence(const double *d_gk, double *outRms)324: void CostFunction::aaConvergence(const double *d_gk, double *outRms)
320: {325: {
321:         using namespace gpu_rigid_bodies;326:         using namespace gpu_rigid_bodies;
322: 327: 
323:         *outRms = 0.0;328:         *outRms = 0.0;
324: 329: 
383:         CudaCheckError();388:         CudaCheckError();
384:         cudaDeviceSynchronize();389:         cudaDeviceSynchronize();
385: 390: 
386:         // Reduction for rigid bodies with more than 32 sites.391:         // Reduction for rigid bodies with more than 32 sites.
387:         if (m_nLargeRigidBody != 0){392:         if (m_nLargeRigidBody != 0){
388:                 int blockSize = 1024;393:                 int blockSize = 1024;
389:                 blockDim.x = blockSize;394:                 blockDim.x = blockSize;
390: 395: 
391:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies.396:                 // Round rigidMaxSite to nearest block size so that per block reduction doesn't sum components from different bodies.
392:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize;397:                 int roundedMaxSite = ((m_rigidMaxSite + blockSize - 1)/blockSize)*blockSize;
 398:                 CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) );
393:                 int numThreads = roundedMaxSite*m_nLargeRigidBody;399:                 int numThreads = roundedMaxSite*m_nLargeRigidBody;
394: 400: 
395:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x;401:                 int blocks = (numThreads + blockDim.x - 1)/blockDim.x;
396:                 gridDim.x = blocks;402:                 gridDim.x = blocks;
397: 403: 
398:                 double *d_outArray;404:                 double3 *d_outArray;
399:                 CudaSafeCall( cudaMalloc(&d_outArray, 3 * blocks * sizeof(double)) );405:                 CudaSafeCall( cudaMalloc(&d_outArray, blocks * sizeof(double)) );
400: 406: 
401:                 gpu_rigid_bodies::fullReduce1b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups,407:                 gpu_rigid_bodies::fullReduce1b<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, m_d_nRigidSitesPerBody, m_d_rigidGroups,
402:                                 m_d_rigidMaxSite, m_d_largeRigidIndices, d_outArray, d_tempArray, roundedMaxSite);408:                                 m_d_rigidMaxSite, m_d_largeRigidIndices, d_outArray, d_tempArray);
403:                 CudaCheckError(); 
404:                 cudaDeviceSynchronize(); 
405: 409: 
406:                 while (blocks > m_nLargeRigidBody) {410:                 while (blocks > m_nLargeRigidBody) {
407:                         int outSize = blocks;411:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::outSize, &blocks, sizeof(int)) );
408: 412: 
409:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize;413:                         roundedMaxSite = ((blocks/m_nLargeRigidBody + blockSize - 1)/blockSize)*blockSize;
 414:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::roundMaxSite, &roundedMaxSite, sizeof(int)) );
410:                         numThreads = roundedMaxSite*m_nLargeRigidBody;415:                         numThreads = roundedMaxSite*m_nLargeRigidBody;
411: 416: 
412:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x;417:                         blocks = (numThreads + blockDim.x - 1)/blockDim.x;
413:                         gridDim.x = blocks;418:                         gridDim.x = blocks;
 419:                         CudaSafeCall( cudaMemcpyToSymbol(gpu_rigid_bodies::nBlocks, &blocks, sizeof(int)) );
414: 420: 
415:                         gpu_rigid_bodies::fullReduce2c<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, d_torques, 421:                         gpu_rigid_bodies::fullReduce2c<<<gridDim, blockDim>>>(m_d_nLargeRigidBody, d_torques, 
416:                                         m_d_largeRigidIndices, d_outArray, m_d_nRigidBody, roundedMaxSite, outSize, blocks);422:                                         m_d_largeRigidIndices, d_outArray, m_d_nRigidBody);
417:                         CudaCheckError(); 
418:                         cudaDeviceSynchronize(); 
419:                 }423:                 }
420: 424: 
421:                 CudaSafeCall( cudaFree(d_outArray) );425:                 CudaSafeCall( cudaFree(d_outArray) );
422:         }426:         }
423: 427: 
424:         CudaSafeCall( cudaFree(d_tempArray) );428:         CudaSafeCall( cudaFree(d_tempArray) );
425: 429: 
426:         double *d_rmsArray;430:         double *d_rmsArray;
427:         CudaSafeCall( cudaMalloc(&d_rmsArray, m_nRigidBody * sizeof(double)) );431:         CudaSafeCall( cudaMalloc(&d_rmsArray, m_nRigidBody * sizeof(double)) );
428: 432: 
441:         double temp2 = 0.0;445:         double temp2 = 0.0;
442: 446: 
443:         // Reduction. 447:         // Reduction. 
444:         blockDim.x = 512;448:         blockDim.x = 512;
445:         gridDim.x = min((m_nRigidBody + blockDim.x - 1) / blockDim.x, 1024);449:         gridDim.x = min((m_nRigidBody + blockDim.x - 1) / blockDim.x, 1024);
446: 450: 
447:         double *d_outArray;451:         double *d_outArray;
448:         CudaSafeCall( cudaMalloc(&d_outArray, gridDim.x * sizeof(double)) );452:         CudaSafeCall( cudaMalloc(&d_outArray, gridDim.x * sizeof(double)) );
449: 453: 
450:         gpu_rigid_bodies::fullReduce<<<gridDim, blockDim>>>(d_rmsArray, d_outArray, m_nRigidBody);454:         gpu_rigid_bodies::fullReduce<<<gridDim, blockDim>>>(d_rmsArray, d_outArray, m_nRigidBody);
451:         CudaCheckError(); 
452:         cudaDeviceSynchronize(); 
453:  
454:         gpu_rigid_bodies::fullReduce<<<1, 1024>>>(d_outArray, d_outArray, gridDim.x);455:         gpu_rigid_bodies::fullReduce<<<1, 1024>>>(d_outArray, d_outArray, gridDim.x);
455:         CudaCheckError();456:         CudaCheckError();
456:         cudaDeviceSynchronize();457:         cudaDeviceSynchronize();
457: 458: 
458:         CudaSafeCall( cudaMemcpy(&temp1, d_outArray, sizeof(double), cudaMemcpyDeviceToHost) );459:         CudaSafeCall( cudaMemcpy(&temp1, d_outArray, sizeof(double), cudaMemcpyDeviceToHost) );
459: 460: 
460:         CudaSafeCall( cudaFree(d_outArray) );461:         CudaSafeCall( cudaFree(d_outArray) );
461:         CudaSafeCall( cudaFree(d_rmsArray) );462:         CudaSafeCall( cudaFree(d_rmsArray) );
462: 463: 
463:         if (m_nDegFreedom > 6*m_nRigidBody) {464:         if (m_nDegFreedom > 6*m_nRigidBody) {
464:                 int arraysize = m_nDegFreedom - 6*m_nRigidBody;465:                 int arraysize = m_nDegFreedom - 6*m_nRigidBody;
465: 466: 
466:                 // Dot product.  467:                 // Dot product.  
467:                 blockDim.x = 512;468:                 blockDim.x = 512;
468:                 gridDim.x = min((arraysize + blockDim.x - 1) / blockDim.x, 1024);469:                 gridDim.x = min((arraysize + blockDim.x - 1) / blockDim.x, 1024);
469: 470: 
470:                 CudaSafeCall( cudaMalloc(&d_outArray, gridDim.x * sizeof(double)) );471:                 CudaSafeCall( cudaMalloc(&d_outArray, gridDim.x * sizeof(double)) );
471: 472: 
472:                 gpu_rigid_bodies::fullDotReduce<<<gridDim, blockDim>>>((m_d_gkRigid + 6*m_nRigidBody), d_outArray, arraysize);473:                 gpu_rigid_bodies::fullDotReduce<<<gridDim, blockDim>>>((m_d_gkRigid + 6*m_nRigidBody), d_outArray, arraysize);
473:                 CudaCheckError(); 
474:                 cudaDeviceSynchronize(); 
475:  
476:                 gpu_rigid_bodies::fullReduce<<<1, 1024>>>(d_outArray, d_outArray, gridDim.x);474:                 gpu_rigid_bodies::fullReduce<<<1, 1024>>>(d_outArray, d_outArray, gridDim.x);
477:                 CudaCheckError();475:                 CudaCheckError();
478:                 cudaDeviceSynchronize();476:                 cudaDeviceSynchronize();
479: 477: 
480:                 CudaSafeCall( cudaMemcpy(&temp2, d_outArray, sizeof(double), cudaMemcpyDeviceToHost) );478:                 CudaSafeCall( cudaMemcpy(&temp2, d_outArray, sizeof(double), cudaMemcpyDeviceToHost) );
481: 479: 
482:                 CudaSafeCall( cudaFree(d_outArray) );480:                 CudaSafeCall( cudaFree(d_outArray) );
483:         }481:         }
484: 482: 
485:         *outRms = temp1 + temp2;483:         *outRms = temp1 + temp2;
783:                                         d_grmi[5+9*thisRigidBody]*781:                                         d_grmi[5+9*thisRigidBody]*
784:                                         m_d_sitesRigidBody[i+(*m_d_rigidMaxSite)+thisRigidBody*3*(*m_d_rigidMaxSite)] + 782:                                         m_d_sitesRigidBody[i+(*m_d_rigidMaxSite)+thisRigidBody*3*(*m_d_rigidMaxSite)] + 
785:                                         d_grmi[8+9*thisRigidBody]*783:                                         d_grmi[8+9*thisRigidBody]*
786:                                         m_d_sitesRigidBody[i+2*(*m_d_rigidMaxSite)+thisRigidBody*3*(*m_d_rigidMaxSite)];784:                                         m_d_sitesRigidBody[i+2*(*m_d_rigidMaxSite)+thisRigidBody*3*(*m_d_rigidMaxSite)];
787:                         }785:                         }
788: 786: 
789:                         tid += blockDim.x * gridDim.x;787:                         tid += blockDim.x * gridDim.x;
790:                 }788:                 }
791:         }789:         }
792: 790: 
793:         __global__ void singleAtoms(double *d_x, const double *m_d_xRigid, const int *m_d_rigidSingles, const int *m_d_nRigidBody, 791:         __global__ void singleAtoms(double *d_x, const double *m_d_xRigid, const int *m_d_rigidSingles, const int *m_d_nRigidBody)
794:                         const int singlesThreads) 
795:         {792:         {
796:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;793:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
797: 794: 
798:                 while (tid < singlesThreads) {795:                 while (tid < singlesThreads) {
799:                         int thisAtom = m_d_rigidSingles[tid];796:                         int thisAtom = m_d_rigidSingles[tid];
800: 797: 
801:                         d_x[3*thisAtom-3] = m_d_xRigid[6*(*m_d_nRigidBody)+3*tid+0];798:                         d_x[3*thisAtom-3] = m_d_xRigid[6*(*m_d_nRigidBody)+3*tid+0];
802:                         d_x[3*thisAtom-2] = m_d_xRigid[6*(*m_d_nRigidBody)+3*tid+1];799:                         d_x[3*thisAtom-2] = m_d_xRigid[6*(*m_d_nRigidBody)+3*tid+1];
803:                         d_x[3*thisAtom-1] = m_d_xRigid[6*(*m_d_nRigidBody)+3*tid+2];800:                         d_x[3*thisAtom-1] = m_d_xRigid[6*(*m_d_nRigidBody)+3*tid+2];
804: 801: 
1001:                                         d_gk[3*myAtom-2]*dr3.y + 998:                                         d_gk[3*myAtom-2]*dr3.y + 
1002:                                         d_gk[3*myAtom-1]*dr3.z;999:                                         d_gk[3*myAtom-1]*dr3.z;
1003: 1000: 
1004:                         }1001:                         }
1005: 1002: 
1006:                         tid += blockDim.x * gridDim.x;1003:                         tid += blockDim.x * gridDim.x;
1007:                 }1004:                 }
1008:         }1005:         }
1009: 1006: 
1010:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles, 1007:         __global__ void gradSingleAtoms(const double *d_gk, double *m_d_gkRigid, const int *m_d_rigidSingles, 
1011:                         const int *m_d_nRigidBody, const int singlesThreads)1008:                         const int *m_d_nRigidBody)
1012:         {1009:         {
1013:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1010:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1014: 1011: 
1015:                 while (tid < singlesThreads) {1012:                 while (tid < singlesThreads) {
1016:                         int thisAtom = m_d_rigidSingles[tid];1013:                         int thisAtom = m_d_rigidSingles[tid];
1017: 1014: 
1018:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+0] = d_gk[3*thisAtom-3];1015:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+0] = d_gk[3*thisAtom-3];
1019:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+1] = d_gk[3*thisAtom-2];1016:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+1] = d_gk[3*thisAtom-2];
1020:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+2] = d_gk[3*thisAtom-1];1017:                         m_d_gkRigid[6*(*m_d_nRigidBody)+3*tid+2] = d_gk[3*thisAtom-1];
1021: 1018: 
1337: 1334: 
1338:                         //Final reduce within first warp1335:                         //Final reduce within first warp
1339:                         if (wid==0) {1336:                         if (wid==0) {
1340:                                 val = warpReduceSum(val);1337:                                 val = warpReduceSum(val);
1341:                         } 1338:                         } 
1342: 1339: 
1343:                         return val;1340:                         return val;
1344:                 }1341:                 }
1345: 1342: 
1346:         __global__ void fullReduce1a(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 1343:         __global__ void fullReduce1a(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
1347:                         const int *m_d_rigidMaxSite, const double *d_gk, const int *m_d_largeRigidIndices, double *d_outArray, 1344:                         const int *m_d_rigidMaxSite, const double *d_gk, const int *m_d_largeRigidIndices, double3 *d_outArray) 
1348:                         const int roundedMaxSite)  
1349:         {1345:         {
1350:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1346:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1351: 1347: 
1352:                 while (tid < (roundedMaxSite*(*m_d_nLargeRigidBody))) {1348:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) {
1353: 1349: 
1354:                         double3 elements;1350:                         double3 elements;
1355:                         elements.x = 0.0;1351:                         elements.x = 0.0;
1356:                         elements.y = 0.0;1352:                         elements.y = 0.0;
1357:                         elements.z = 0.0;1353:                         elements.z = 0.0;
1358: 1354: 
1359:                         int index = ((tid + roundedMaxSite) % roundedMaxSite);1355:                         int index = ((tid + roundMaxSite) % roundMaxSite);
1360:                         int thisBody = tid/roundedMaxSite;1356:                         int thisBody = tid/roundMaxSite;
1361:                         int currentBody = m_d_largeRigidIndices[thisBody];1357:                         int currentBody = m_d_largeRigidIndices[thisBody];
1362: 1358: 
1363:                         if (index < m_d_nRigidSitesPerBody[currentBody]) {1359:                         if (index < m_d_nRigidSitesPerBody[currentBody]) {
1364:                                 int myAtom = m_d_rigidGroups[index+(*m_d_rigidMaxSite)*currentBody];1360:                                 int myAtom = m_d_rigidGroups[index+(*m_d_rigidMaxSite)*currentBody];
1365:                                 elements.x = d_gk[3*myAtom-3];1361:                                 elements.x = d_gk[3*myAtom-3];
1366:                                 elements.y = d_gk[3*myAtom-2];1362:                                 elements.y = d_gk[3*myAtom-2];
1367:                                 elements.z = d_gk[3*myAtom-1];1363:                                 elements.z = d_gk[3*myAtom-1];
1368:                         }1364:                         }
1369: 1365: 
1370:                         __syncthreads(); 
1371:  
1372:                         elements.x = blockReduceSum(elements.x);1366:                         elements.x = blockReduceSum(elements.x);
1373:                         elements.y = blockReduceSum(elements.y);1367:                         elements.y = blockReduceSum(elements.y);
1374:                         elements.z = blockReduceSum(elements.z);1368:                         elements.z = blockReduceSum(elements.z);
1375: 1369: 
1376:                         __syncthreads(); 
1377:  
1378:                         if (threadIdx.x==0) {1370:                         if (threadIdx.x==0) {
1379:                                 d_outArray[3*blockIdx.x+0] = elements.x;1371:                                 d_outArray[blockIdx.x].x = elements.x;
1380:                                 d_outArray[3*blockIdx.x+1] = elements.y;1372:                                 d_outArray[blockIdx.x].y = elements.y;
1381:                                 d_outArray[3*blockIdx.x+2] = elements.z;1373:                                 d_outArray[blockIdx.x].z = elements.z;
1382:                         }1374:                         }
1383: 1375: 
1384:                         tid += blockDim.x * gridDim.x;1376:                         tid += blockDim.x * gridDim.x;
1385:                 }1377:                 }
1386:         }1378:         }
1387: 1379: 
1388:         __global__ void fullReduce1b(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 1380:         __global__ void fullReduce1b(const int *m_d_nLargeRigidBody, const int *m_d_nRigidSitesPerBody, const int *m_d_rigidGroups, 
1389:                         const int *m_d_rigidMaxSite, const int *m_d_largeRigidIndices, double *d_outArray, 1381:                         const int *m_d_rigidMaxSite, const int *m_d_largeRigidIndices, double3 *d_outArray, 
1390:                         const double *d_tempArray, const int roundedMaxSite) 1382:                         const double *d_tempArray) 
1391:         {1383:         {
1392:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1384:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1393: 1385: 
1394:                 while (tid < (roundedMaxSite*(*m_d_nLargeRigidBody))) {1386:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) {
1395: 1387: 
1396:                         double3 elements;1388:                         double3 elements;
1397:                         elements.x = 0.0;1389:                         elements.x = 0.0;
1398:                         elements.y = 0.0;1390:                         elements.y = 0.0;
1399:                         elements.z = 0.0;1391:                         elements.z = 0.0;
1400: 1392: 
1401:                         int index = ((tid + roundedMaxSite) % roundedMaxSite);1393:                         int index = ((tid + roundMaxSite) % roundMaxSite);
1402:                         int thisBody = tid/roundedMaxSite;1394:                         int thisBody = tid/roundMaxSite;
1403:                         int currentBody = m_d_largeRigidIndices[thisBody];1395:                         int currentBody = m_d_largeRigidIndices[thisBody];
1404: 1396: 
1405:                         if (index < m_d_nRigidSitesPerBody[currentBody]) {1397:                         if (index < m_d_nRigidSitesPerBody[currentBody]) {
1406:                                 elements.x = d_tempArray[0+3*index+3*(*m_d_rigidMaxSite)*currentBody];1398:                                 elements.x = d_tempArray[0+3*index+3*(*m_d_rigidMaxSite)*currentBody];
1407:                                 elements.y = d_tempArray[1+3*index+3*(*m_d_rigidMaxSite)*currentBody];1399:                                 elements.y = d_tempArray[1+3*index+3*(*m_d_rigidMaxSite)*currentBody];
1408:                                 elements.z = d_tempArray[2+3*index+3*(*m_d_rigidMaxSite)*currentBody];1400:                                 elements.z = d_tempArray[2+3*index+3*(*m_d_rigidMaxSite)*currentBody];
1409:                         }1401:                         }
1410: 1402: 
1411:                         __syncthreads(); 
1412:  
1413:                         elements.x = blockReduceSum(elements.x);1403:                         elements.x = blockReduceSum(elements.x);
1414:                         elements.y = blockReduceSum(elements.y);1404:                         elements.y = blockReduceSum(elements.y);
1415:                         elements.z = blockReduceSum(elements.z);1405:                         elements.z = blockReduceSum(elements.z);
1416: 1406: 
1417:                         __syncthreads(); 
1418:  
1419:                         if (threadIdx.x==0) {1407:                         if (threadIdx.x==0) {
1420:                                 d_outArray[3*blockIdx.x+0] = elements.x;1408:                                 d_outArray[blockIdx.x].x = elements.x;
1421:                                 d_outArray[3*blockIdx.x+1] = elements.y;1409:                                 d_outArray[blockIdx.x].y = elements.y;
1422:                                 d_outArray[3*blockIdx.x+2] = elements.z;1410:                                 d_outArray[blockIdx.x].z = elements.z;
1423:                         }1411:                         }
1424: 1412: 
1425:                         tid += blockDim.x * gridDim.x;1413:                         tid += blockDim.x * gridDim.x;
1426:                 }1414:                 }
1427:         }1415:         }
1428: 1416: 
1429:         __global__ void fullReduce2a(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, 1417:         __global__ void fullReduce2a(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, double3 *d_outArray) 
1430:                         double *d_outArray, const int roundedMaxSite, const int outSize, const int blocks)  
1431:         {1418:         {
1432:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1419:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1433: 1420: 
1434:                 while (tid < (roundedMaxSite*(*m_d_nLargeRigidBody))) {1421:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) {
1435: 1422: 
1436:                         double3 elements;1423:                         double3 elements;
1437:                         elements.x = 0.0;1424:                         elements.x = 0.0;
1438:                         elements.y = 0.0;1425:                         elements.y = 0.0;
1439:                         elements.z = 0.0;1426:                         elements.z = 0.0;
1440: 1427: 
1441:                         int index = ((tid + roundedMaxSite) % roundedMaxSite);1428:                         int index = ((tid + roundMaxSite) % roundMaxSite);
1442:                         int thisBody = tid/roundedMaxSite;1429:                         int thisBody = tid/roundMaxSite;
1443:                         int currentBody = m_d_largeRigidIndices[thisBody];1430:                         int currentBody = m_d_largeRigidIndices[thisBody];
1444: 1431: 
1445:                         int sectionSize = outSize/(*m_d_nLargeRigidBody);1432:                         int sectionSize = outSize/(*m_d_nLargeRigidBody);
1446:                         if (index < sectionSize) {1433:                         if (index < sectionSize) {
1447:                                 int var = index+thisBody*sectionSize;1434:                                 elements.x = d_outArray[index+thisBody*sectionSize].x;
1448:                                 elements.x = d_outArray[3*var+0];1435:                                 elements.y = d_outArray[index+thisBody*sectionSize].y;
1449:                                 elements.y = d_outArray[3*var+1];1436:                                 elements.z = d_outArray[index+thisBody*sectionSize].z;
1450:                                 elements.z = d_outArray[3*var+2]; 
1451:                         }1437:                         }
1452: 1438: 
1453:                         __syncthreads(); 
1454:  
1455:                         elements.x = blockReduceSum(elements.x);1439:                         elements.x = blockReduceSum(elements.x);
1456:                         elements.y = blockReduceSum(elements.y);1440:                         elements.y = blockReduceSum(elements.y);
1457:                         elements.z = blockReduceSum(elements.z);1441:                         elements.z = blockReduceSum(elements.z);
1458: 1442: 
1459:                         __syncthreads(); 
1460:  
1461:                         if (threadIdx.x==0) {1443:                         if (threadIdx.x==0) {
1462:                                 if (blocks == (*m_d_nLargeRigidBody)) {1444:                                 if (nBlocks == (*m_d_nLargeRigidBody)) {
1463:                                         m_d_gkRigid[3*currentBody] = elements.x;1445:                                         m_d_gkRigid[3*currentBody] = elements.x;
1464:                                         m_d_gkRigid[3*currentBody+1] = elements.y;1446:                                         m_d_gkRigid[3*currentBody+1] = elements.y;
1465:                                         m_d_gkRigid[3*currentBody+2] = elements.z;1447:                                         m_d_gkRigid[3*currentBody+2] = elements.z;
1466:                                 }1448:                                 }
1467:                                 else {1449:                                 else {
1468:                                         d_outArray[3*blockIdx.x+0] = elements.x;1450:                                         d_outArray[blockIdx.x].x = elements.x;
1469:                                         d_outArray[3*blockIdx.x+1] = elements.y;1451:                                         d_outArray[blockIdx.x].y = elements.y;
1470:                                         d_outArray[3*blockIdx.x+2] = elements.z;1452:                                         d_outArray[blockIdx.x].z = elements.z;
1471:                                 }1453:                                 }
1472:                         }1454:                         }
1473: 1455: 
1474:                         tid += blockDim.x * gridDim.x;1456:                         tid += blockDim.x * gridDim.x;
1475:                 }1457:                 }
1476:         }1458:         }
1477: 1459: 
1478:         __global__ void fullReduce2b(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, 1460:         __global__ void fullReduce2b(const int *m_d_nLargeRigidBody, double *m_d_gkRigid, const int *m_d_largeRigidIndices, 
1479:                         double *d_outArray, const int *m_d_nRigidBody, const int roundedMaxSite, const int outSize, const int blocks) 1461:                         double3 *d_outArray, const int *m_d_nRigidBody) 
1480:         {1462:         {
1481:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1463:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1482: 1464: 
1483:                 while (tid < (roundedMaxSite*(*m_d_nLargeRigidBody))) {1465:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) {
1484: 1466: 
1485:                         double3 elements;1467:                         double3 elements;
1486:                         elements.x = 0.0;1468:                         elements.x = 0.0;
1487:                         elements.y = 0.0;1469:                         elements.y = 0.0;
1488:                         elements.z = 0.0;1470:                         elements.z = 0.0;
1489: 1471: 
1490:                         int index = ((tid + roundedMaxSite) % roundedMaxSite);1472:                         int index = ((tid + roundMaxSite) % roundMaxSite);
1491:                         int thisBody = tid/roundedMaxSite;1473:                         int thisBody = tid/roundMaxSite;
1492:                         int currentBody = m_d_largeRigidIndices[thisBody];1474:                         int currentBody = m_d_largeRigidIndices[thisBody];
1493: 1475: 
1494:                         int sectionSize = outSize/(*m_d_nLargeRigidBody);1476:                         int sectionSize = outSize/(*m_d_nLargeRigidBody);
1495:                         if (index < sectionSize) {1477:                         if (index < sectionSize) {
1496:                                 int var = index+thisBody*sectionSize;1478:                                 elements.x = d_outArray[index+thisBody*sectionSize].x;
1497:                                 elements.x = d_outArray[3*var+0];1479:                                 elements.y = d_outArray[index+thisBody*sectionSize].y;
1498:                                 elements.y = d_outArray[3*var+1];1480:                                 elements.z = d_outArray[index+thisBody*sectionSize].z;
1499:                                 elements.z = d_outArray[3*var+2]; 
1500:                         }1481:                         }
1501: 1482: 
1502:                         __syncthreads(); 
1503:  
1504:                         elements.x = blockReduceSum(elements.x);1483:                         elements.x = blockReduceSum(elements.x);
1505:                         elements.y = blockReduceSum(elements.y);1484:                         elements.y = blockReduceSum(elements.y);
1506:                         elements.z = blockReduceSum(elements.z);1485:                         elements.z = blockReduceSum(elements.z);
1507: 1486: 
1508:                         __syncthreads(); 
1509:  
1510:                         if (threadIdx.x==0) {1487:                         if (threadIdx.x==0) {
1511:                                 if (blocks == (*m_d_nLargeRigidBody)) {1488:                                 if (nBlocks == (*m_d_nLargeRigidBody)) {
1512:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody] = elements.x;1489:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody] = elements.x;
1513:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+1] = elements.y;1490:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+1] = elements.y;
1514:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+2] = elements.z;1491:                                         m_d_gkRigid[3*(*m_d_nRigidBody)+3*currentBody+2] = elements.z;
1515:                                 }1492:                                 }
1516:                                 else {1493:                                 else {
1517:                                         d_outArray[3*blockIdx.x+0] = elements.x;1494:                                         d_outArray[blockIdx.x].x = elements.x;
1518:                                         d_outArray[3*blockIdx.x+1] = elements.y;1495:                                         d_outArray[blockIdx.x].y = elements.y;
1519:                                         d_outArray[3*blockIdx.x+2] = elements.z;1496:                                         d_outArray[blockIdx.x].z = elements.z;
1520:                                 }1497:                                 }
1521:                         }1498:                         }
1522: 1499: 
1523:                         tid += blockDim.x * gridDim.x;1500:                         tid += blockDim.x * gridDim.x;
1524:                 }1501:                 }
1525:         }1502:         }
1526: 1503: 
1527:         __global__ void fullReduce2c(const int *m_d_nLargeRigidBody, double *d_torques, const int *m_d_largeRigidIndices, 1504:         __global__ void fullReduce2c(const int *m_d_nLargeRigidBody, double *d_torques, const int *m_d_largeRigidIndices, double3 *d_outArray, const int *m_d_nRigidBody) 
1528:                         double *d_outArray, const int *m_d_nRigidBody, const int roundedMaxSite, const int outSize,  
1529:                         const int blocks)  
1530:         {1505:         {
1531:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;1506:                 int tid = blockIdx.x * blockDim.x + threadIdx.x;
1532: 1507: 
1533:                 while (tid < (roundedMaxSite*(*m_d_nLargeRigidBody))) {1508:                 while (tid < (roundMaxSite*(*m_d_nLargeRigidBody))) {
1534: 1509: 
1535:                         double3 elements;1510:                         double3 elements;
1536:                         elements.x = 0.0;1511:                         elements.x = 0.0;
1537:                         elements.y = 0.0;1512:                         elements.y = 0.0;
1538:                         elements.z = 0.0;1513:                         elements.z = 0.0;
1539: 1514: 
1540:                         int index = ((tid + roundedMaxSite) % roundedMaxSite); // Technically unnecessary with one body at a time - equivalent to tid1515:                         int index = ((tid + roundMaxSite) % roundMaxSite); // Technically unnecessary with one body at a time - equivalent to tid
1541:                         int thisBody = tid/roundedMaxSite; // Also unnecessary with one body - could just pass as argument1516:                         int thisBody = tid/roundMaxSite; // Also unnecessary with one body - could just pass as argument
1542:                         int currentBody = m_d_largeRigidIndices[thisBody];1517:                         int currentBody = m_d_largeRigidIndices[thisBody];
1543: 1518: 
1544:                         int sectionSize = outSize/(*m_d_nLargeRigidBody);1519:                         int sectionSize = outSize/(*m_d_nLargeRigidBody);
1545:                         if (index < sectionSize) {1520:                         if (index < sectionSize) {
1546:                                 int var = index+thisBody*sectionSize;1521:                                 elements.x = d_outArray[index+thisBody*sectionSize].x;
1547:                                 elements.x = d_outArray[3*var+0];1522:                                 elements.y = d_outArray[index+thisBody*sectionSize].y;
1548:                                 elements.y = d_outArray[3*var+1];1523:                                 elements.z = d_outArray[index+thisBody*sectionSize].z;
1549:                                 elements.z = d_outArray[3*var+2]; 
1550:                         }1524:                         }
1551: 1525: 
1552:                         __syncthreads(); 
1553:  
1554:                         elements.x = blockReduceSum(elements.x);1526:                         elements.x = blockReduceSum(elements.x);
1555:                         elements.y = blockReduceSum(elements.y);1527:                         elements.y = blockReduceSum(elements.y);
1556:                         elements.z = blockReduceSum(elements.z);1528:                         elements.z = blockReduceSum(elements.z);
1557: 1529: 
1558:                         __syncthreads(); 
1559:  
1560:                         if (threadIdx.x==0) {1530:                         if (threadIdx.x==0) {
1561:                                 if (blocks == (*m_d_nLargeRigidBody)) {1531:                                 if (nBlocks == (*m_d_nLargeRigidBody)) {
1562:                                         d_torques[3*currentBody] = elements.x;1532:                                         d_torques[3*currentBody] = elements.x;
1563:                                         d_torques[3*currentBody+1] = elements.y;1533:                                         d_torques[3*currentBody+1] = elements.y;
1564:                                         d_torques[3*currentBody+2] = elements.z;1534:                                         d_torques[3*currentBody+2] = elements.z;
1565:                                 }1535:                                 }
1566:                                 else {1536:                                 else {
1567:                                         d_outArray[3*blockIdx.x+0] = elements.x;1537:                                         d_outArray[blockIdx.x].x = elements.x;
1568:                                         d_outArray[3*blockIdx.x+1] = elements.y;1538:                                         d_outArray[blockIdx.x].y = elements.y;
1569:                                         d_outArray[3*blockIdx.x+2] = elements.z;1539:                                         d_outArray[blockIdx.x].z = elements.z;
1570:                                 }1540:                                 }
1571:                         }1541:                         }
1572: 1542: 
1573:                         tid += blockDim.x * gridDim.x;1543:                         tid += blockDim.x * gridDim.x;
1574:                 }1544:                 }
1575: 1545: 
1576:         }1546:         }
1577: 1547: 
1578:         __global__ void fullReduce(const double *in, double* out, const int N) 1548:         __global__ void fullReduce(double *in, double* out, int N) 
1579:         {1549:         {
1580:                 double sum = 0;1550:                 double sum = 0;
1581:                 //reduce multiple elements per thread1551:                 //reduce multiple elements per thread
1582:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) {1552:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) {
1583:                         sum += in[i];1553:                         sum += in[i];
1584:                 }1554:                 }
1585:                 sum = blockReduceSum(sum);1555:                 sum = blockReduceSum(sum);
1586:                 if (threadIdx.x==0) {1556:                 if (threadIdx.x==0) {
1587:                         out[blockIdx.x]=sum;1557:                         out[blockIdx.x]=sum;
1588:                 }1558:                 }
1589:         }1559:         }
1590: 1560: 
1591:         __global__ void fullDotReduce(const double *in, double* out, const int N)1561:         __global__ void fullDotReduce(double *in, double* out, int N)
1592:         {1562:         {
1593:                 double sum = 0;1563:                 double sum = 0;
1594:                 //reduce multiple elements per thread1564:                 //reduce multiple elements per thread
1595:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) {1565:                 for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) {
1596:                         sum += in[i]*in[i];1566:                         sum += in[i]*in[i];
1597:                 }1567:                 }
1598:                 sum = blockReduceSum(sum);1568:                 sum = blockReduceSum(sum);
1599:                 if (threadIdx.x==0) {1569:                 if (threadIdx.x==0) {
1600:                         out[blockIdx.x]=sum;1570:                         out[blockIdx.x]=sum;
1601:                 }1571:                 }


r32161/setup_bfgsts.cu 2017-03-22 12:30:27.764236529 +0000 r32160/setup_bfgsts.cu 2017-03-22 12:30:31.296283126 +0000
111: 111: 
112: extern "C" void setup_bfgsts(int *nAtoms, double *coords, double *xGradEps, _Bool *isConverged, double *energy, int *xMaxIter, 112: extern "C" void setup_bfgsts(int *nAtoms, double *coords, double *xGradEps, _Bool *isConverged, double *energy, int *xMaxIter, 
113:                 int *bItMax, int *bIter, double *xMaxStep, double *xMaxFkRise, double *outRms, char *cudaPotential, 113:                 int *bItMax, int *bIter, double *xMaxStep, double *xMaxFkRise, double *outRms, char *cudaPotential, 
114:                 _Bool *printDebug, _Bool *timeCuda, int *nPotCalls, _Bool *isColdFusion, double *coldFusionLim, double *xH0, 114:                 _Bool *printDebug, _Bool *timeCuda, int *nPotCalls, _Bool *isColdFusion, double *coldFusionLim, double *xH0, 
115:                 int *xUpdates, double *evalMin, double *evec, _Bool *isAtomisticNotRigid, int *nDegFreedomF, int *nRigidBodyF, 115:                 int *xUpdates, double *evalMin, double *evec, _Bool *isAtomisticNotRigid, int *nDegFreedomF, int *nRigidBodyF, 
116:                 int *nRigidSitesPerBodyF, int *rigidGroupsF, int *rigidMaxSiteF, double *sitesRigidBodyF, int *rigidSinglesF, 116:                 int *nRigidSitesPerBodyF, int *rigidGroupsF, int *rigidMaxSiteF, double *sitesRigidBodyF, int *rigidSinglesF, 
117:                 double *rigidInverseF, int *nSecDiag, double *evalPercentEps, double *pushOffCut, double *pushOffMag, 117:                 double *rigidInverseF, int *nSecDiag, double *evalPercentEps, double *pushOffCut, double *pushOffMag, 
118:                 double *maxEvecStep, double *maxMaxStep, double *minMaxStep, double *trustRadius, int *pMaxIter1, int *pMaxIter2, 118:                 double *maxEvecStep, double *maxMaxStep, double *minMaxStep, double *trustRadius, int *pMaxIter1, int *pMaxIter2, 
119:                 double *evecOverlapTol, int *pUpdates, double *gradEps, double *pMaxStep, double *pMaxFkRise, double *pH0, 119:                 double *evecOverlapTol, int *pUpdates, double *gradEps, double *pMaxStep, double *pMaxFkRise, double *pH0, 
120:                 double *evStepMagTol, _Bool *shouldFreeze, _Bool *isAtomFrozen, int *nFreeze, double *potTimeElapsed, 120:                 double *evStepMagTol, _Bool *shouldFreeze, _Bool *isAtomFrozen, int *nFreeze, double *potTimeElapsed, 
121:                 bool *isAaConvergence, int *nAddTarget, double *ljAddRep, double *ljAddAtt)121:                 bool *isAaConvergence)
122: {122: {
123:         const size_t numDimensions = 3*(*nAtoms);123:         const size_t numDimensions = 3*(*nAtoms);
124:         const double aaConvThreshold = 5*(*gradEps);124:         const double aaConvThreshold = 5*(*gradEps);
125: 125: 
126:         // Rigid body parameters may not be intialised in Fortran code if rigid body framework not being used. 126:         // Rigid body parameters may not be intialised in Fortran code if rigid body framework not being used. 
127:         int nDegFreedom  = 0;127:         int nDegFreedom  = 0;
128:         int nRigidBody   = 0;128:         int nRigidBody   = 0;
129:         int rigidMaxSite = 0;129:         int rigidMaxSite = 0;
130: 130: 
131:         int *nRigidSitesPerBody = NULL;131:         int *nRigidSitesPerBody = NULL;
151:         if (*printDebug) {151:         if (*printDebug) {
152:                 debugPrinting.setPrintingOn();152:                 debugPrinting.setPrintingOn();
153:         }153:         }
154: 154: 
155:         Timer timer_potential     ("GPU_potential");155:         Timer timer_potential     ("GPU_potential");
156:         timer_potential.setTimingOn();156:         timer_potential.setTimingOn();
157: 157: 
158:         // Set up cuBLAS. 158:         // Set up cuBLAS. 
159:         Cublas cublas;159:         Cublas cublas;
160: 160: 
161:         // L specifies the Lennard-Jones potential. 161:         // L specifies the Lennard-Jones potential for up to 1024 atoms. 
162:         if (*cudaPotential == 'L') {162:         if (*cudaPotential == 'L') {
 163:                 const size_t atomMax = 1024;
 164:                 if (numDimensions > 3*atomMax) {
 165:                         std::cerr << "Lennard-Jones is currently only supported for a maximum of 1024 atoms. " << std::endl;
 166:                         exit(EXIT_FAILURE);
 167:                 }
163:                 // Create an instance of the appropriate class for the potential, LjPotential. 168:                 // Create an instance of the appropriate class for the potential, LjPotential. 
164:                 LjPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 169:                 LjPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 
165:                                 nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, coords, 170:                                 nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, coords, 
166:                                 *nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, isAtomFrozen, 171:                                 *nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, isAtomFrozen, 
167:                                 *nFreeze, *isAaConvergence, *nAddTarget, ljAddRep, ljAddAtt);172:                                 *nFreeze, *isAaConvergence);
168: 173: 
169:                 setup<LjPotential>(debugPrinting, cublas, numDimensions, xGradEps, isConverged, energy, xMaxIter, bIter, xMaxStep, xMaxFkRise, 174:                 setup<LjPotential>(debugPrinting, cublas, numDimensions, xGradEps, isConverged, energy, xMaxIter, bIter, xMaxStep, xMaxFkRise, 
170:                                 outRms, potential, printDebug, timeCuda, nPotCalls, isColdFusion, coldFusionLim, xH0, xUpdates, 175:                                 outRms, potential, printDebug, timeCuda, nPotCalls, isColdFusion, coldFusionLim, xH0, xUpdates, 
171:                                 evec, evalPercentEps, evalMin, coords, pushOffCut, pushOffMag, maxEvecStep, maxMaxStep, 176:                                 evec, evalPercentEps, evalMin, coords, pushOffCut, pushOffMag, maxEvecStep, maxMaxStep, 
172:                                 minMaxStep, trustRadius, pMaxIter1, pMaxIter2, bItMax, evecOverlapTol, pUpdates, gradEps, 177:                                 minMaxStep, trustRadius, pMaxIter1, pMaxIter2, bItMax, evecOverlapTol, pUpdates, gradEps, 
173:                                 pMaxStep, pMaxFkRise, pH0, evStepMagTol);178:                                 pMaxStep, pMaxFkRise, pH0, evStepMagTol);
174:         }179:         }
175:         // A specifies the AMBER potential. 180:         // A specifies the AMBER potential. 
176:         else if (*cudaPotential == 'A') {181:         else if (*cudaPotential == 'A') {
177:                 AmberPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 182:                 AmberPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, rigidMaxSite, 
178:                                 nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, coords, 183:                                 nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, coords, 
179:                                 *nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, isAtomFrozen, 184:                                 *nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, isAtomFrozen, 
180:                                 *nFreeze, *isAaConvergence, *nAddTarget, ljAddRep, ljAddAtt);185:                                 *nFreeze, *isAaConvergence);
181: 186: 
182:                 setup<AmberPotential>(debugPrinting, cublas, numDimensions, xGradEps, isConverged, energy, xMaxIter, bIter, xMaxStep, xMaxFkRise, 187:                 setup<AmberPotential>(debugPrinting, cublas, numDimensions, xGradEps, isConverged, energy, xMaxIter, bIter, xMaxStep, xMaxFkRise, 
183:                                 outRms, potential, printDebug, timeCuda, nPotCalls, isColdFusion, coldFusionLim, xH0, 188:                                 outRms, potential, printDebug, timeCuda, nPotCalls, isColdFusion, coldFusionLim, xH0, 
184:                                 xUpdates, evec, evalPercentEps, evalMin, coords, pushOffCut, pushOffMag, maxEvecStep, 189:                                 xUpdates, evec, evalPercentEps, evalMin, coords, pushOffCut, pushOffMag, maxEvecStep, 
185:                                 maxMaxStep, minMaxStep, trustRadius, pMaxIter1, pMaxIter2, bItMax, evecOverlapTol, pUpdates, 190:                                 maxMaxStep, minMaxStep, trustRadius, pMaxIter1, pMaxIter2, bItMax, evecOverlapTol, pUpdates, 
186:                                 gradEps, pMaxStep, pMaxFkRise, pH0, evStepMagTol);191:                                 gradEps, pMaxStep, pMaxFkRise, pH0, evStepMagTol);
187:         }192:         }
188:         else {193:         else {
189:                 std::cerr << "The specified potential has not been recognised. " << std::endl;194:                 std::cerr << "The specified potential has not been recognised. " << std::endl;
190:                 exit(EXIT_FAILURE);195:                 exit(EXIT_FAILURE);


r32161/setup_lbfgs.cu 2017-03-22 12:30:28.864251041 +0000 r32160/setup_lbfgs.cu 2017-03-22 12:30:32.400297692 +0000
 74: } 74: }
 75:  75: 
 76:  76: 
 77:  77: 
 78: extern "C" void setup_lbfgs(int *n, double *x, double *gradEps, _Bool *isConverged, double *energy, int *maxIter, int *itDone,  78: extern "C" void setup_lbfgs(int *n, double *x, double *gradEps, _Bool *isConverged, double *energy, int *maxIter, int *itDone, 
 79:                 double *maxStep, double *maxFkRise, double *outRms, char *cudaPot, _Bool *printDebug, _Bool *timeCuda,  79:                 double *maxStep, double *maxFkRise, double *outRms, char *cudaPot, _Bool *printDebug, _Bool *timeCuda, 
 80:                 int *nPotCalls, _Bool *isColdFusion, double *coldFusionLim, double *H0, int *updates,  80:                 int *nPotCalls, _Bool *isColdFusion, double *coldFusionLim, double *H0, int *updates, 
 81:                 _Bool *isAtomisticNotRigid, int *nDegFreedomF, int *nRigidBodyF, int *nRigidSitesPerBodyF,  81:                 _Bool *isAtomisticNotRigid, int *nDegFreedomF, int *nRigidBodyF, int *nRigidSitesPerBodyF, 
 82:                 int *rigidGroupsF, int *rigidMaxSiteF, double *sitesRigidBodyF, int *rigidSinglesF,  82:                 int *rigidGroupsF, int *rigidMaxSiteF, double *sitesRigidBodyF, int *rigidSinglesF, 
 83:                 double *aaConvThresholdF, double *rigidInverseF, _Bool *projectGrad, _Bool *shouldFreeze,  83:                 double *aaConvThresholdF, double *rigidInverseF, _Bool *projectGrad, _Bool *shouldFreeze, 
 84:                 _Bool *isAtomFrozen, int *nFreeze, double *potTimeElapsed, bool *isAaConvergence, int *nAddTarget, double *ljAddRep,  84:                 _Bool *isAtomFrozen, int *nFreeze, double *potTimeElapsed, bool *isAaConvergence)
 85:                 double *ljAddAtt) 
 86: { 85: {
 87:         const size_t numDimensions = *n; 86:         const size_t numDimensions = *n;
 88:  87: 
 89:         const int nSecDiag = 0; 88:         const int nSecDiag = 0;
 90:  89: 
 91:         const double aaConvThreshold = 5*(*aaConvThresholdF); 90:         const double aaConvThreshold = 5*(*aaConvThresholdF);
 92:  91: 
 93:         // Rigid body parameters may not be intialised if rigid body framework not being used. 92:         // Rigid body parameters may not be intialised if rigid body framework not being used.
 94:         int nDegFreedom  = 0; 93:         int nDegFreedom  = 0;
 95:         int nRigidBody   = 0; 94:         int nRigidBody   = 0;
120:         if (*printDebug) {119:         if (*printDebug) {
121:                 debugPrinting.setPrintingOn();120:                 debugPrinting.setPrintingOn();
122:         }121:         }
123: 122: 
124:         Timer timer_potential     ("GPU_potential");123:         Timer timer_potential     ("GPU_potential");
125:         timer_potential.setTimingOn();124:         timer_potential.setTimingOn();
126: 125: 
127:         // Set up cuBLAS. 126:         // Set up cuBLAS. 
128:         Cublas cublas;127:         Cublas cublas;
129: 128: 
130:         // L specifies the Lennard-Jones potential.129:         // L specifies the Lennard-Jones potential for up to 1024 atoms.
131:         if (*cudaPot == 'L') {130:         if (*cudaPot == 'L') {
 131:                 const size_t atom_max = 1024;
 132:                 if (numDimensions > 3*atom_max) {
 133:                         std::cerr << "Lennard-Jones is currently only supported for a maximum of 1024 atoms. " << std::endl;
 134:                         exit(EXIT_FAILURE);
 135:                 }
132:                 // Create an instance of the appropriate class for the potential, LjPotential.136:                 // Create an instance of the appropriate class for the potential, LjPotential.
133:                 LjPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, 137:                 LjPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, 
134:                                 rigidMaxSite, nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, 138:                                 rigidMaxSite, nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, 
135:                                 coords, nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, 139:                                 coords, nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, 
136:                                 isAtomFrozen, *nFreeze, *isAaConvergence, *nAddTarget, ljAddRep, ljAddAtt);140:                                 isAtomFrozen, *nFreeze, *isAaConvergence);
137: 141: 
138:                 setup<LjPotential>(debugPrinting, cublas, numDimensions, x, gradEps, isConverged, energy, maxIter, itDone, 142:                 setup<LjPotential>(debugPrinting, cublas, numDimensions, x, gradEps, isConverged, energy, maxIter, itDone, 
139:                                 maxStep, maxFkRise, outRms, potential, printDebug, timeCuda, nPotCalls, H0, updates, isColdFusion, 143:                                 maxStep, maxFkRise, outRms, potential, printDebug, timeCuda, nPotCalls, H0, updates, isColdFusion, 
140:                                 projectGrad);144:                                 projectGrad);
141:         }145:         }
142:         // A specifies the AMBER potential.146:         // A specifies the AMBER potential.
143:         else if (*cudaPot == 'A') {147:         else if (*cudaPot == 'A') {
144:                 AmberPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, 148:                 AmberPotential potential(debugPrinting, timer_potential, cublas, numDimensions, nDegFreedom, nRigidBody, 
145:                                 rigidMaxSite, nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, 149:                                 rigidMaxSite, nRigidSitesPerBody, rigidGroups, sitesRigidBody, rigidSingles, rigidInverse, 
146:                                 coords, nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, 150:                                 coords, nSecDiag, *isAtomisticNotRigid, aaConvThreshold, *coldFusionLim, *shouldFreeze, 
147:                                 isAtomFrozen, *nFreeze, *isAaConvergence, *nAddTarget, ljAddRep, ljAddAtt);151:                                 isAtomFrozen, *nFreeze, *isAaConvergence);
148: 152: 
149:                 setup<AmberPotential>(debugPrinting, cublas, numDimensions, x, gradEps, isConverged, energy, maxIter, itDone, 153:                 setup<AmberPotential>(debugPrinting, cublas, numDimensions, x, gradEps, isConverged, energy, maxIter, itDone, 
150:                                 maxStep, maxFkRise, outRms, potential, printDebug, timeCuda, nPotCalls, H0, updates, isColdFusion, 154:                                 maxStep, maxFkRise, outRms, potential, printDebug, timeCuda, nPotCalls, H0, updates, isColdFusion, 
151:                                 projectGrad); 155:                                 projectGrad); 
152:         }156:         }
153:         else {157:         else {
154:                 std::cerr << "The specified potential has not been recognised" << std::endl;158:                 std::cerr << "The specified potential has not been recognised" << std::endl;
155:                 exit(EXIT_FAILURE);159:                 exit(EXIT_FAILURE);
156:         }160:         }
157: 161: 


r32161/sparse_bench.F90 2017-03-22 12:30:30.408271411 +0000 r32160/sparse_bench.F90 2017-03-22 12:30:33.960318274 +0000
  1: SUBROUTINE BENCH_SPARSE(COORDINATES, QUENCH_NUM, HESS_CUTOFF, POT_ENERGY, LOG_PROD)  1: SUBROUTINE BENCH_SPARSE(COORDINATES, QUENCH_NUM, HESS_CUTOFF, POT_ENERGY, LOG_PROD)
  2:    USE COMMONS, ONLY: NATOMS  2:    USE COMMONS, ONLY: NATOMS, DEBUG
  3:    USE MODHESS  3:    USE MODHESS
  4: #ifdef __SPARSE  4: #ifdef __SPARSE
  5:    USE MODSPARSEHESS  5:    USE MODSPARSEHESS
  6:    USE SHIFT_HESS  6:    USE SHIFT_HESS
  7: #endif  7: #endif
  8:    IMPLICIT NONE  8:    IMPLICIT NONE
  9: ! Arguments  9: ! Arguments
 10:    DOUBLE PRECISION, INTENT(IN)  :: COORDINATES(3*NATOMS) 10:    DOUBLE PRECISION, INTENT(IN)  :: COORDINATES(3*NATOMS)
 11:    INTEGER, INTENT(IN)           :: QUENCH_NUM 11:    INTEGER, INTENT(IN)           :: QUENCH_NUM
 12:    DOUBLE PRECISION, INTENT(IN)  :: HESS_CUTOFF 12:    DOUBLE PRECISION, INTENT(IN)  :: HESS_CUTOFF
 13:    DOUBLE PRECISION, INTENT(OUT) :: POT_ENERGY 13:    DOUBLE PRECISION, INTENT(OUT) :: POT_ENERGY
 14:    DOUBLE PRECISION, INTENT(OUT) :: LOG_PROD 14:    DOUBLE PRECISION, INTENT(OUT) :: LOG_PROD
 15: ! Timing variables 15: ! Timing variables
 16:    DOUBLE PRECISION              :: TIME_DSYEV, TIME_FILT, TIME_SPARSE 16:    DOUBLE PRECISION              :: TIME_DSYEV, TIME_FILT, TIME_SPARSE
 17:    DOUBLE PRECISION              :: T_START, T_END 17:    DOUBLE PRECISION              :: T_START, T_END
 18: ! File variables 18: ! File variables
 19:    INTEGER                       :: SUMMARY_FILE 19:    INTEGER                       :: SUMMARY_FILE, EFILE1, EFILE2, EFILE3
 20:    INTEGER                       :: GETUNIT 20:    INTEGER                       :: GETUNIT
 21: ! DSYEV variables 21: ! DSYEV variables
 22:    INTEGER, PARAMETER            :: LWORK = 100000 22:    INTEGER, PARAMETER            :: LWORK = 100000
 23:    DOUBLE PRECISION              :: WORK(LWORK) 23:    DOUBLE PRECISION              :: WORK(LWORK)
 24:    INTEGER                       :: INFO 24:    INTEGER                       :: INFO
 25: ! Other variables 25: ! Other variables
 26:    DOUBLE PRECISION              :: LOG_PROD_DSYEV, LOG_PROD_FILT, LOG_PROD_SPARSE 26:    DOUBLE PRECISION              :: LOG_PROD_DSYEV, LOG_PROD_FILT, LOG_PROD_SPARSE
 27:    DOUBLE PRECISION              :: EVALUES(3*NATOMS) 27:    DOUBLE PRECISION              :: EVALUES(3*NATOMS)
 28:    DOUBLE PRECISION              :: SHIFT_ARRAY(6) 28:    DOUBLE PRECISION              :: SHIFT_ARRAY(6)
 29:    DOUBLE PRECISION              :: GRAD(3*NATOMS) 29:    DOUBLE PRECISION              :: GRAD(3*NATOMS)
 30:    DOUBLE PRECISION, ALLOCATABLE :: HESS_COPY(:, :) 30:    DOUBLE PRECISION, ALLOCATABLE :: HESS_COPY(:, :), COORDSCOPY(:)
 31:    INTEGER                       :: NUM_ZERO_EVS  31:    INTEGER                       :: NUM_ZERO_EVS , MYVAR, I, J
  32: 
  33:    INTEGER, PARAMETER              :: HESS_OUT1 = 4739
 32:  34: 
 33: #ifdef __SPARSE 35: #ifdef __SPARSE
 34: ! Assume 6 zero eigenvalues 36: ! Assume 6 zero eigenvalues
 35:    NUM_ZERO_EVS = 6 37:    NUM_ZERO_EVS = 6
 36:  38: 
 37: ! Calculate hessian and mass weight 39: ! Calculate hessian and mass weight
 38:    IF (ALLOCATED(HESS)) DEALLOCATE(HESS) 40:    IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
 39:    ALLOCATE(HESS(3*NATOMS, 3*NATOMS)) 41:    ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
 40:    CALL POTENTIAL(COORDINATES, GRAD, POT_ENERGY, .TRUE., .TRUE.) 42:    CALL POTENTIAL(COORDINATES, GRAD, POT_ENERGY, .TRUE., .TRUE.)
 41:    CALL MASSWT() 43:    CALL MASSWT()
 42:  44: 
 43: ! Copy the original (unfiltered) hessian to somewhere safe 45: ! Copy the original (unfiltered) hessian to somewhere safe
 44:    IF (ALLOCATED(HESS_COPY)) DEALLOCATE(HESS_COPY) 46:    IF (ALLOCATED(HESS_COPY)) DEALLOCATE(HESS_COPY)
 45:    ALLOCATE(HESS_COPY(3*NATOMS, 3*NATOMS)) 47:    ALLOCATE(HESS_COPY(3*NATOMS, 3*NATOMS))
 46:    HESS_COPY = HESS 48:    HESS_COPY = HESS
 47:  49: 
  50: ! Copy the original coordinates to somewhere safe
  51:    IF (ALLOCATED(COORDSCOPY)) DEALLOCATE(COORDSCOPY)
  52:    ALLOCATE(COORDSCOPY(3*NATOMS))
  53:    COORDSCOPY = COORDINATES
  54: 
 48: ! Assign shift array 55: ! Assign shift array
 49:    SHIFT_ARRAY(1:6) = 1.0D0 56:    SHIFT_ARRAY(1:6) = 1.0D0
 50:  57: 
 51: ! Sparse section 58: ! Sparse section
 52:    CALL CPU_TIME(T_START) 59:    CALL CPU_TIME(T_START)
 53:    CALL SHIFT_HESS_ZEROS(COORDINATES, SHIFT_ARRAY) 60:    CALL SHIFT_HESS_ZEROS(COORDINATES, SHIFT_ARRAY)
 54:    CALL FILTER_ZEROS(HESS, HESS_CUTOFF) 61:    CALL FILTER_ZEROS(HESS, HESS_CUTOFF)
 55:    CALL GET_DETERMINANT(3*NATOMS, LOG_PROD_SPARSE, QUENCH_NUM) 62:    CALL GET_DETERMINANT(3*NATOMS, LOG_PROD_SPARSE, QUENCH_NUM)
 56:    CALL CPU_TIME(T_END) 63:    CALL CPU_TIME(T_END)
 57:    TIME_SPARSE = T_END - T_START 64:    TIME_SPARSE = T_END - T_START
 65:    TIME_DSYEV = T_END - T_START 72:    TIME_DSYEV = T_END - T_START
 66:    LOG_PROD_DSYEV = SUM(DLOG(EVALUES(NUM_ZERO_EVS + 1:))) 73:    LOG_PROD_DSYEV = SUM(DLOG(EVALUES(NUM_ZERO_EVS + 1:)))
 67:  74: 
 68: ! Need to pass back the log product of frequencies for actual FEBH (sqrt of eigenvalues) 75: ! Need to pass back the log product of frequencies for actual FEBH (sqrt of eigenvalues)
 69:    LOG_PROD = 0.5D0 * LOG_PROD_DSYEV 76:    LOG_PROD = 0.5D0 * LOG_PROD_DSYEV
 70:  77: 
 71: ! Print summary 78: ! Print summary
 72:    SUMMARY_FILE = GETUNIT()  79:    SUMMARY_FILE = GETUNIT() 
 73:    OPEN(UNIT=SUMMARY_FILE, FILE='sparse_benchmarks', STATUS='UNKNOWN', ACCESS='APPEND') 80:    OPEN(UNIT=SUMMARY_FILE, FILE='sparse_benchmarks', STATUS='UNKNOWN', ACCESS='APPEND')
 74:    WRITE(SUMMARY_FILE, '(A, I12)') 'Quench', QUENCH_NUM 81:    WRITE(SUMMARY_FILE, '(A, I12)') 'Quench', QUENCH_NUM
 75:    WRITE(SUMMARY_FILE, '(A, A, E18.10, A, F10.4, A, E18.10, A, F8.2)') 'DSYEV             ||', & 82:    WRITE(SUMMARY_FILE, '(A, A, E18.10, A, F10.4, A, E18.10, A, F8.2)') 'DSYEV             !!||', &
 76:                                                                     ' log_prod= ', LOG_PROD_DSYEV, & 83:                                                                     ' log_prod= ', LOG_PROD_DSYEV, &
 77:                                                                     ' time= ', TIME_DSYEV, & 84:                                                                     ' time= ', TIME_DSYEV, &
 78:                                                                     ' error= ', LOG_PROD_DSYEV - LOG_PROD_DSYEV, & 85:                                                                     ' error= ', LOG_PROD_DSYEV - LOG_PROD_DSYEV, &
 79:                                                                     ' % error= ', 100.0D0 * (LOG_PROD_DSYEV - LOG_PROD_DSYEV) / LOG_PROD_DSYEV 86:                                                                     ' % error= ', 100.0D0 * (LOG_PROD_DSYEV - LOG_PROD_DSYEV) / LOG_PROD_DSYEV
 80:    WRITE(SUMMARY_FILE, '(A, A, E18.10, A, F10.4, A, E18.10, A, F8.2)') 'Sparse            ||', & 87:    WRITE(SUMMARY_FILE, '(A, A, E18.10, A, F10.4, A, E18.10, A, F8.2)') 'Sparse            ||', &
 81:                                                                     ' log_prod= ', LOG_PROD_SPARSE, & 88:                                                                     ' log_prod= ', LOG_PROD_SPARSE, &
 82:                                                                     ' time= ', TIME_SPARSE, & 89:                                                                     ' time= ', TIME_SPARSE, &
 83:                                                                     ' error= ', LOG_PROD_SPARSE - LOG_PROD_DSYEV, & 90:                                                                     ' error= ', LOG_PROD_SPARSE - LOG_PROD_DSYEV, &
 84:                                                                     ' % error= ', 100.0D0 * (LOG_PROD_SPARSE - LOG_PROD_DSYEV) / LOG_PROD_DSYEV 91:                                                                     ' % error= ', 100.0D0 * (LOG_PROD_SPARSE - LOG_PROD_DSYEV) / LOG_PROD_DSYEV
 85:    CLOSE(SUMMARY_FILE)  92:    CLOSE(SUMMARY_FILE) 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0