Added diagnostic method to GBVI

platforms/reference/src/ReferenceKernels.cpp

@@ -859,11 +859,14 @@ ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
 }
 
 void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii ) {
+
     int numParticles = system.getNumParticles();
+
     charges.resize(numParticles);
     vector<RealOpenMM> atomicRadii(numParticles);
     vector<RealOpenMM> scaledRadii(numParticles);
     vector<RealOpenMM> gammas(numParticles);
+
     for (int i = 0; i < numParticles; ++i) {
         double charge, radius, gamma;
         force.getParticleParameters(i, charge, radius, gamma);
@@ -872,15 +875,19 @@ void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIFo
         gammas[i]        = static_cast<RealOpenMM>(gamma);
         scaledRadii[i]   = static_cast<RealOpenMM>(inputScaledRadii[i]);
     }
+
     GBVIParameters * gBVIParameters = new GBVIParameters(numParticles);
+
     gBVIParameters->setAtomicRadii(atomicRadii);
     gBVIParameters->setGammaParameters(gammas);
     gBVIParameters->setScaledRadii(scaledRadii);
     gBVIParameters->setSolventDielectric(static_cast<RealOpenMM>(force.getSolventDielectric()));
     gBVIParameters->setSoluteDielectric(static_cast<RealOpenMM>(force.getSoluteDielectric()));
+
     gBVIParameters->setBornRadiusScalingMethod(force.getBornRadiusScalingMethod());
     gBVIParameters->setQuinticUpperBornRadiusLimit(static_cast<RealOpenMM>(force.getQuinticUpperBornRadiusLimit()));
     gBVIParameters->setQuinticLowerLimitFactor(static_cast<RealOpenMM>(force.getQuinticLowerLimitFactor()));
+
     if (force.getNonbondedMethod() != GBVIForce::NoCutoff)
         gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
     isPeriodic = (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic);

platforms/reference/src/gbsa/CpuGBVI.cpp

@@ -724,7 +724,7 @@ void CpuGBVI::computeBornForces( std::vector<RealVec>& atomCoordinates, const Re
         }
     }
 
-    //printGbvi( atomCoordinates, partialCharges, bornRadii, bornForces, forces, "Post loop2", stderr );
+    //printGbvi( atomCoordinates, partialCharges, bornRadii, bornForces, forces, "GBVI: Post loop2", stderr );
 
     // convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
 
@@ -776,26 +776,48 @@ void CpuGBVI::printGbvi( const std::vector<OpenMM::RealVec>& atomCoordinates, co
     const RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
     RealOpenMM tau                            = static_cast<RealOpenMM>(gbviParameters->getTau());  
 
+    int useComparisonFormat                   = 1;
+
     (void) fprintf( log, "Reference Gbvi     %s atoms=%d\n", idString.c_str(), numberOfAtoms );
     (void) fprintf( log, "    tau            %15.7e\n", tau ); 
     (void) fprintf( log, "    scaleMethod    %d (QuinticEnum=%d)\n", 
                     _gbviParameters->getBornRadiusScalingMethod(), GBVIParameters::QuinticSpline );
     (void) fprintf( log, "    preFactor      %15.7e)\n", preFactor );
  
-    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-        (void) fprintf( log, "%6d r=%15.7e rSc=%15.7e swd=%15.7e tau*gam=%15.7e q=%15.7e", atomI,
-                        atomicRadii[atomI],    
-                        scaledRadii[atomI],    
-                        switchDeriviative[atomI],    
-                        tau*gammaParameters[atomI],    
-                        partialCharges[atomI] );
-        if( bornRadii.size() > atomI ){
-             (void) fprintf( log, " bR=%15.7e", bornRadii[atomI] );
-        }
-        if( bornForces.size() > atomI ){
-             (void) fprintf( log, " tau*bF=%15.7e", tau*bornForces[atomI] );    
-        }
-        (void) fprintf( log, "\n" );
+    if( useComparisonFormat ){
+        (void) fprintf( log, "  br bF swd r scR tau*gamma q)\n" );
+        for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+            (void) fprintf( log, "%6d ", atomI );
+            if( bornRadii.size() > atomI ){
+                 (void) fprintf( log, "%15.7e ", bornRadii[atomI] );
+            }
+            if( bornForces.size() > atomI ){
+                 (void) fprintf( log, "%15.7e ", tau*bornForces[atomI] );    
+            }
+            (void) fprintf( log, " %15.7e %15.7e %15.7e %15.7e %15.7e",
+                            switchDeriviative[atomI],    
+                            atomicRadii[atomI],    
+                            scaledRadii[atomI],    
+                            tau*gammaParameters[atomI],    
+                            partialCharges[atomI] );
+            (void) fprintf( log, "\n" );
+        }   
+    } else {
+        for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+            (void) fprintf( log, "%6d r=%15.7e rSc=%15.7e swd=%15.7e tau*gam=%15.7e q=%15.7e", atomI,
+                            atomicRadii[atomI],    
+                            scaledRadii[atomI],    
+                            switchDeriviative[atomI],    
+                            tau*gammaParameters[atomI],    
+                            partialCharges[atomI] );
+            if( bornRadii.size() > atomI ){
+                 (void) fprintf( log, " bR=%15.7e", bornRadii[atomI] );
+            }
+            if( bornForces.size() > atomI ){
+                 (void) fprintf( log, " tau*bF=%15.7e", tau*bornForces[atomI] );    
+            }
+            (void) fprintf( log, "\n" );
+        }   
     }   
 
     return;