Converted units from kcal/A to Joule/nm

platforms/reference/src/gbsa/CpuImplicitSolvent.cpp

@@ -699,8 +699,7 @@ int CpuImplicitSolvent::computeImplicitSolventForces( RealOpenMM** atomCoordinat
 
    // compute forces
 
-   computeBornEnergyForces( getBornRadii(), atomCoordinates,
-                                                partialCharges, forces );
+   computeBornEnergyForces( getBornRadii(), atomCoordinates, partialCharges, forces );
 
    // diagnostics
 

platforms/reference/src/gbsa/CpuObc.cpp

@@ -225,7 +225,7 @@ int CpuObc::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadi
    static const RealOpenMM half    = (RealOpenMM) 0.5;
    static const RealOpenMM fourth  = (RealOpenMM) 0.25;
 
-   // static const char* methodName = "\nCpuObc::computeBornRadii";
+   static const char* methodName = "\nCpuObc::computeBornRadii";
 
    // ---------------------------------------------------------------------------------------
 
@@ -320,8 +320,7 @@ if( logFile && atomI == 0 ){
 if( logFile && atomI >= 0 ){
    (void) fprintf( logFile, "\nRRQ %d sum=%12.6e tanhS=%12.6e radI=%.5f %.5f born=%12.6e obc=%12.6e",
                    atomI, sum, tanhSum, radiusI, offsetRadiusI, bornRadii[atomI], obcChain[atomI] );
-} 
-*/
+} */
 
    }
 
@@ -381,7 +380,6 @@ int CpuObc::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoo
    // constants
 
    const RealOpenMM preFactor           = obcParameters->getPreFactor();
-   const RealOpenMM electricConstant    = obcParameters->getElectricConstant();
    const RealOpenMM dielectricOffset    = obcParameters->getDielectricOffset();
 
    // ---------------------------------------------------------------------------------------
@@ -594,7 +592,15 @@ int CpuObc::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoo
       obcChainTemp[atomI]        = (one - tanhSum*tanhSum)*obcChainTemp[atomI]/radiusI;
    }
 
-   setEnergy( obcEnergy );
+   // cal to Joule conversion
+
+   RealOpenMM conversion = 0.4184;  
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+      forces[atomI][0] *= conversion;
+      forces[atomI][1] *= conversion;
+      forces[atomI][2] *= conversion;
+   }
+   setEnergy( obcEnergy*conversion );
 
    // copy new Born radii and obcChain values into permanent array
 
@@ -882,9 +888,9 @@ int CpuObc::computeBornEnergyForcesPrint( RealOpenMM* bornRadii, RealOpenMM** at
 #pragma warning(disable:4996)
 #endif
 
-//FILE* logFile = NULL;
+FILE* logFile = NULL;
 //FILE* logFile = SimTKOpenMMLog::getSimTKOpenMMLogFile( );
-FILE* logFile = fopen( "bF", "w" );
+//FILE* logFile = fopen( "bF", "w" );
 
 #if defined(_MSC_VER)
 #pragma warning(pop)
@@ -895,7 +901,6 @@ FILE* logFile = fopen( "bF", "w" );
    // constants
 
    const RealOpenMM preFactor           = obcParameters->getPreFactor();
-   const RealOpenMM electricConstant    = obcParameters->getElectricConstant();
    const RealOpenMM dielectricOffset    = obcParameters->getDielectricOffset();
 
    // ---------------------------------------------------------------------------------------
@@ -919,6 +924,14 @@ FILE* logFile = fopen( "bF", "w" );
     
    if( includeAceApproximation() ){
       computeAceNonPolarForce( obcParameters, bornRadii, &obcEnergy, bornForces );
+
+      if( logFile ){
+         (void) fprintf( logFile, "\nACE E=%.5e\n", obcEnergy );
+         for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+            (void) fprintf( logFile, "   %d bR=%.6e bF=%.6e\n", atomI, bornRadii[atomI], bornForces[atomI] );
+         }
+      }
+
    }
 
    // ---------------------------------------------------------------------------------------
@@ -989,7 +1002,9 @@ FILE* logFile = fopen( "bF", "w" );
          obcEnergy         += Gpol;
          bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
 
-if( logFile && (atomI == -1 || atomJ == -1) ){
+/*
+if( logFile ){
+//if( logFile && (atomI == -1 || atomJ == -1) ){
 //   (void) fprintf( logFile, "\nWWX %d %d F[%.6e %.6e %.6e] bF=[%.6e %.6e] Gpl[%.6e %.6e %.6e] rb[%6.4f %7.4f] rs[%6.4f %7.4f] ",
 //                    atomI, atomJ,
 //                    forces[atomI][0],  forces[atomI][1],  forces[atomI][2],
@@ -1004,16 +1019,17 @@ if( logFile && (atomI == -1 || atomJ == -1) ){
                     bornRadii[atomI], bornRadii[atomJ],
                     dGpol_dalpha2_ij*bornRadii[atomJ], dGpol_dalpha2_ij*bornRadii[atomI] );
 }
+*/
       }
 
 
    }
 
 if( logFile ){
-   (void) fprintf( logFile, "\nWXX bF & F E=%.8e", obcEnergy );
+   (void) fprintf( logFile, "\nWXX bF & F E=%.8e preFactor=%.5f", obcEnergy, preFactor );
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-      (void) fprintf( logFile, "\nWXX %d %.6e q=%.3f F[%.6e %.6e %.6e] ",
-                      atomI, partialCharges[atomI],  bornForces[atomI], forces[atomI][0],  forces[atomI][1],  forces[atomI][2] );
+      (void) fprintf( logFile, "\nWXX %d q=%.4f bR=%.5e bF=%.3f F[%.6e %.6e %.6e] ",
+                      atomI, partialCharges[atomI], bornRadii[atomI],  bornForces[atomI], forces[atomI][0],  forces[atomI][1],  forces[atomI][2] );
    }
 }
 
@@ -1187,7 +1203,13 @@ if( logFile && atomI >= 0 ){
 
    }
 
-   setEnergy( obcEnergy );
+   RealOpenMM conversion = 0.4184;  
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+      forces[atomI][0] *= conversion;
+      forces[atomI][1] *= conversion;
+      forces[atomI][2] *= conversion;
+   }
+   setEnergy( obcEnergy*conversion );
 
    if( 0 ){
 

platforms/reference/src/gbsa/ImplicitSolventParameters.cpp

@@ -478,14 +478,12 @@ int ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
    Set AtomicRadii array
 
    @param atomicRadii array of atomic radii
-   @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                  SimTKOpenMMCommon::MdUnits 
 
    @return SimTKOpenMMCommon::DefaultReturn
 
    --------------------------------------------------------------------------------------- */
 
-int ImplicitSolventParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii, int units ){
+int ImplicitSolventParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -511,17 +509,8 @@ int ImplicitSolventParameters::setAtomicRadii( const RealOpenMMVector& atomicRad
       numberOfAtoms = numberOfAtoms < (int) atomicRadii.size() ? numberOfAtoms : (int) atomicRadii.size();
    }
 
-   // force kcal/A units
-
-   if( units == SimTKOpenMMCommon::MdUnits ){
-      RealOpenMM ten = (RealOpenMM) 10.0;
-      for( int ii = 0; ii < numberOfAtoms; ii++ ){
-         _atomicRadii[ii] = ten*atomicRadii[ii];
-      }
-   } else {
-      for( int ii = 0; ii < numberOfAtoms; ii++ ){
-         _atomicRadii[ii] = atomicRadii[ii];
-      }
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      _atomicRadii[ii] = atomicRadii[ii];
    }
 
    return SimTKOpenMMCommon::DefaultReturn;
@@ -532,14 +521,12 @@ int ImplicitSolventParameters::setAtomicRadii( const RealOpenMMVector& atomicRad
    Set AtomicRadii array
 
    @param atomicRadii array of atomic radii
-   @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                  SimTKOpenMMCommon::MdUnits 
 
    @return SimTKOpenMMCommon::DefaultReturn
 
    --------------------------------------------------------------------------------------- */
-
-int ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
+/*
+int ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -559,19 +546,12 @@ int ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii, int unit
    
    // force kcal/A units
 
-   if( units == SimTKOpenMMCommon::MdUnits ){
-      RealOpenMM ten = (RealOpenMM) 10.0;
-      for( int ii = 0; ii < numberOfAtoms; ii++ ){
-         _atomicRadii[ii] = ten*atomicRadii[ii];
-      }
-   } else {
-      for( int ii = 0; ii < numberOfAtoms; ii++ ){
-         _atomicRadii[ii] = atomicRadii[ii];
-      }
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      _atomicRadii[ii] = atomicRadii[ii];
    }
 
    return SimTKOpenMMCommon::DefaultReturn;
-}
+} */
 
 /**---------------------------------------------------------------------------------------
 
@@ -656,13 +636,16 @@ void ImplicitSolventParameters::_initializeImplicitSolventConstants( void ){
    _soluteDielectric        = (RealOpenMM)    1.0;
    _solventDielectric       = (RealOpenMM)   78.3;
    _kcalA_To_kJNm           = (RealOpenMM)    0.4184;
-   _probeRadius             = (RealOpenMM)    1.4;
+   _probeRadius             = (RealOpenMM)    0.14;
    _electricConstant        = (RealOpenMM) -166.02691;
 
-   // _pi4Asolv                = (RealOpenMM) PI_M*4.0*0.0049*1000.0;
+   //_pi4Asolv                = (RealOpenMM) PI_M*4.0*0.0049*1000.0;
    //_pi4Asolv                = (RealOpenMM) PI_M*19.6;
-   // _pi4Asolv                = (RealOpenMM) PI_M*4.0*0.0054;
-   _pi4Asolv                = (RealOpenMM) (PI_M*0.0216);
+   //_pi4Asolv                = (RealOpenMM) PI_M*4.0*0.0054;
+   _pi4Asolv                = (RealOpenMM) (PI_M*0.0216*1000.0);
+   //_pi4Asolv                = (RealOpenMM) -400.71504079;
+   //_pi4Asolv                = (RealOpenMM) 0.0;
+
    
    _resetPreFactor();
 }
@@ -771,7 +754,7 @@ int ImplicitSolventParameters::isNotReady( void ) const {
                                              &stdDev, &minValue, &minIndex,
                                              &maxValue, &maxIndex );
 
-   if( average < 1.0 || average > 10.0 || minValue < 0.5 ){
+   if( average < 0.1 || average > 1.0 || minValue < 0.05 ){
       errors++;
       message << "\n   atomic radii appear not to be set correctly -- radii should be in Angstroms";
       message << "\n   average radius=" << average << " min radius=" << minValue << " at atom index=" << minIndex;
@@ -822,6 +805,7 @@ int ImplicitSolventParameters::isNotReady( void ) const {
       warning++;
       message << "\n   Warning: probe radius=" << getProbeRadius() << " is large.";
    }
+/*
 
    if( getPi4Asolv() <= 0.0 ){
       errors++;
@@ -832,6 +816,7 @@ int ImplicitSolventParameters::isNotReady( void ) const {
       warning++;
       message << "\n   Warning: Pi4Asolv=" << getPi4Asolv() << " is large.";
    }
+*/
 
    if( errors || warning ){
       message << std::endl;

platforms/reference/src/gbsa/ImplicitSolventParameters.h

@@ -310,29 +310,24 @@ class ImplicitSolventParameters {
          Set AtomicRadii array
       
          @param atomicRadii vector of atomic radii
-         @param units       units flag SimTKOpenMMCommon::MdUnits or
-                                       SimTKOpenMMCommon::KcalAngUnits
       
          @return SimTKOpenMMCommon::DefaultReturn
       
          --------------------------------------------------------------------------------------- */
 
-      virtual int setAtomicRadii( const RealOpenMMVector& atomicRadii,
-                                  int units = SimTKOpenMMCommon::MdUnits );
+      virtual int setAtomicRadii( const RealOpenMMVector& atomicRadii );
       
       /**---------------------------------------------------------------------------------------
       
          Set AtomicRadii array
       
          @param atomicRadii array of atomic radii
-         @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                        SimTKOpenMMCommon::MdUnits 
       
          @return SimTKOpenMMCommon::DefaultReturn
       
          --------------------------------------------------------------------------------------- */
       
-      virtual int setAtomicRadii( RealOpenMM* atomicRadii, int units );
+      //virtual int setAtomicRadii( RealOpenMM* atomicRadii );
       
       /**---------------------------------------------------------------------------------------
       
@@ -354,7 +349,6 @@ class ImplicitSolventParameters {
          Print state to log file (Simbios)
          
          @param title               title (optional)
-         @param log                 print state to log file
             
          @return SimTKOpenMMCommon::DefaultReturn
             

platforms/reference/src/gbsa/ObcParameters.cpp

@@ -110,7 +110,7 @@ ObcParameters::ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType
    // ---------------------------------------------------------------------------------------
 
    _obcType                      = obcType;
-   _dielectricOffset             = 0.09f;
+   _dielectricOffset             = 0.009f;
    _ownScaledRadiusFactors       = 0;
    _scaledRadiusFactors          = NULL;
 
@@ -315,14 +315,12 @@ RealOpenMM* ObcParameters::getAtomicRadii( void ) const {
    Set AtomicRadii array
 
    @param atomicRadii array of atomic radii
-   @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                  SimTKOpenMMCommon::MdUnits 
 
    @return SimTKOpenMMCommon::DefaultReturn
 
    --------------------------------------------------------------------------------------- */
 
-int ObcParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
+int ObcParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -330,7 +328,7 @@ int ObcParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
 
    // ---------------------------------------------------------------------------------------
 
-   return ImplicitSolventParameters::setAtomicRadii( atomicRadii, units );
+   return ImplicitSolventParameters::setAtomicRadii( atomicRadii );
 }
 
 /**---------------------------------------------------------------------------------------
@@ -338,14 +336,12 @@ int ObcParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
    Set AtomicRadii array
 
    @param atomicRadii vector of atomic radii
-   @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                  SimTKOpenMMCommon::MdUnits 
 
    @return SimTKOpenMMCommon::DefaultReturn
 
    --------------------------------------------------------------------------------------- */
 
-int ObcParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii, int units ){
+int ObcParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -353,7 +349,7 @@ int ObcParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii, int unit
 
    // ---------------------------------------------------------------------------------------
 
-   return ImplicitSolventParameters::setAtomicRadii( atomicRadii, units );
+   return ImplicitSolventParameters::setAtomicRadii( atomicRadii );
 }
 
 /**---------------------------------------------------------------------------------------
@@ -646,7 +642,7 @@ int ObcParameters::isNotReady( void ) const {
 
    if( average < 0.3 || average > 2.0 || minValue < 0.1 ){
       errors++;
-      message << "\n   scale factors for atomic radii appear not to be set correctly -- radii should be in Angstroms";
+      message << "\n   scale factors for atomic radii appear not to be set correctly -- radii should be in nm";
       message << "\n   average radius=" << average << " min radius=" << minValue << " at atom index=" << minIndex;
    }   
 

platforms/reference/src/gbsa/ObcParameters.h

@@ -232,28 +232,24 @@ class ObcParameters : public ImplicitSolventParameters {
          Set AtomicRadii array
       
          @param atomicRadii array of atomic radii
-         @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                        SimTKOpenMMCommon::MdUnits 
       
          @return SimTKOpenMMCommon::DefaultReturn
       
          --------------------------------------------------------------------------------------- */
 
-      int setAtomicRadii( RealOpenMM* atomicRadii, int units = SimTKOpenMMCommon::MdUnits  );
+      int setAtomicRadii( RealOpenMM* atomicRadii );
 
       /**---------------------------------------------------------------------------------------
       
          Set AtomicRadii array
       
          @param atomicRadii vector of atomic radii
-         @param units       units flag: SimTKOpenMMCommon::KcalAngUnits or
-                                        SimTKOpenMMCommon::MdUnits 
       
          @return SimTKOpenMMCommon::DefaultReturn
       
          --------------------------------------------------------------------------------------- */
 
-      int setAtomicRadii( const RealOpenMMVector& atomicRadii, int units = SimTKOpenMMCommon::MdUnits );
+      int setAtomicRadii( const RealOpenMMVector& atomicRadii );
 
       /**---------------------------------------------------------------------------------------
       

platforms/reference/src/gbsa/cpuObcInterface.cpp

@@ -86,7 +86,7 @@ cpuSetObcParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* obc
 
    ObcParameters* obcParameters  = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
    obcParameters->setScaledRadiusFactors( obcScaleFactors );
-   obcParameters->setAtomicRadii( atomicRadii, SimTKOpenMMCommon::KcalAngUnits );
+   obcParameters->setAtomicRadii( atomicRadii );
 
    // dielectric constants