ComputeGridForce.C

Go to the documentation of this file.
00001 
00007 #include "ComputeGridForce.h"
00008 #include "GridForceGrid.h"
00009 #include "Node.h"
00010 
00011 #define MIN_DEBUG_LEVEL 2
00012 //#define DEBUGM
00013 #include "Debug.h"
00014 
00015 #include "GridForceGrid.inl"
00016 #include "MGridforceParams.h"
00017 
00018 //#define GF_FORCE_OUTPUT
00019 //#define GF_FORCE_OUTPUT_FREQ 100
00020 #define GF_OVERLAPCHECK_FREQ 1000
00021 
00022 
00023 ComputeGridForce::ComputeGridForce(ComputeID c, PatchID pid)
00024     : ComputeHomePatch(c,pid)
00025 {
00026 
00027     reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00028 
00029 }
00030 /*                      END OF FUNCTION ComputeGridForce        */
00031 
00032 
00033 ComputeGridForce::~ComputeGridForce()
00034 {
00035     delete reduction;
00036 }
00037 /*                      END OF FUNCTION ~ComputeGridForce       */
00038 
00039 template <class T> void ComputeGridForce::do_calc(T *grid, int gridnum, FullAtom *p, int numAtoms, Molecule *mol, Force *forces, BigReal &energy, Force &extForce, Tensor &extVirial)
00040 {
00041     Real scale;                 // Scaling factor
00042     Charge charge;              // Charge
00043     Vector dV;
00044     float V;
00045     
00046     Vector gfScale = grid->get_scale();
00047 
00048     //  Loop through and check each atom
00049     for (int i = 0; i < numAtoms; i++) {
00050         if (mol->is_atom_gridforced(p[i].id, gridnum))
00051         {
00052             DebugM(1, "Atom " << p[i].id << " is gridforced\n" << endi);
00053             
00054             mol->get_gridfrc_params(scale, charge, p[i].id, gridnum);
00055             
00056             // Wrap coordinates using grid center
00057             Position pos = grid->wrap_position(p[i].position, homePatch->lattice);
00058             DebugM(1, "pos = " << pos << "\n" << endi);
00059             
00060             // Here's where the action happens
00061             int err = grid->compute_VdV(pos, V, dV);
00062             
00063             if (err) {
00064                 DebugM(2, "V = 0\n" << endi);
00065                 DebugM(2, "dV = 0 0 0\n" << endi);
00066                 continue;  // This means the current atom is outside the potential
00067             }
00068             
00069             //Force force = scale * Tensor::diagonal(gfScale) * (-charge * dV);
00070             Force force = -charge * scale * Vector(gfScale.x * dV.x, gfScale.y * dV.y, gfScale.z * dV.z);
00071             
00072 #ifdef DEBUGM
00073             DebugM(2, "scale = " << scale << " gfScale = " << gfScale << " charge = " << charge << "\n" << endi);
00074             
00075             DebugM(2, "V = " << V << "\n" << endi);
00076             DebugM(2, "dV = " << dV << "\n" << endi);
00077             DebugM(2, "grid = " << gridnum << " force = " << force << " pos = " << pos << " V = " << V << " dV = " << dV << " step = " << homePatch->flags.step << " index = " << p[i].id << "\n" << endi);
00078             
00079             DebugM(1, "transform = " << (int)p[i].transform.i << " "
00080                    << (int)p[i].transform.j << " " << (int)p[i].transform.k << "\n" << endi);
00081             
00082             if (V != V) {
00083                 iout << iWARN << "V is NaN!\natomid = " << p[i].id << " loc = " << p[i].position << " V = " << V << "\n" << endi;
00084             }
00085 #endif
00086             
00087             forces[i] += force;
00088             extForce += force;
00089             Position vpos = homePatch->lattice.reverse_transform(p[i].position, p[i].transform);
00090             
00091             //energy -= force * (vpos - homePatch->lattice.origin());
00092             if (gfScale.x == gfScale.y && gfScale.x == gfScale.z)
00093             {
00094                 // only makes sense when scaling is isotropic
00095                 energy += scale * gfScale.x * (charge * V);
00096                 
00097                 // add something when we're off the grid? I'm thinking no
00098             }
00099             extVirial += outer(force,vpos);
00100         }
00101     }
00102 }
00103 
00104 void ComputeGridForce::doForce(FullAtom* p, Results* r)
00105 {
00106     SimParameters *simParams = Node::Object()->simParameters;
00107     Molecule *mol = Node::Object()->molecule;
00108     
00109     Force *forces = r->f[Results::normal];
00110     BigReal energy = 0;
00111     Force extForce = 0.;
00112     Tensor extVirial;
00113     
00114     int numAtoms = homePatch->getNumAtoms();
00115 
00116     if ( mol->numGridforceGrids < 1 ) NAMD_bug("No grids loaded in ComputeGridForce::doForce()");
00117     
00118     for (int gridnum = 0; gridnum < mol->numGridforceGrids; gridnum++) {
00119         GridforceGrid *grid = mol->get_gridfrc_grid(gridnum);
00120         
00121         if (homePatch->flags.step % GF_OVERLAPCHECK_FREQ == 0) {
00122             // only check on node 0 and every GF_OVERLAPCHECK_FREQ steps
00123           if (simParams->langevinPistonOn || simParams->berendsenPressureOn) {
00124                 // check for grid overlap if pressure control is on
00125                 // not needed without pressure control, since the check is also performed on startup
00126       if (!grid->fits_lattice(homePatch->lattice)) {
00127         char errmsg[512];
00128         if (grid->get_checksize()) {
00129           sprintf(errmsg, "Warning: Periodic cell basis too small for Gridforce grid %d.  Set gridforcechecksize off in configuration file to ignore.\n", gridnum);
00130           NAMD_die(errmsg);      
00131         }
00132       }
00133          }
00134         }
00135         
00136         Position center = grid->get_center();
00137         
00138         if (homePatch->flags.step % 100 == 1) {
00139             DebugM(3, "center = " << center << "\n" << endi);
00140             DebugM(3, "e = " << grid->get_e() << "\n" << endi);
00141         }
00142         
00143         if (grid->get_grid_type() == GridforceGrid::GridforceGridTypeFull) {
00144             GridforceFullMainGrid *g = (GridforceFullMainGrid *)grid;
00145             do_calc(g, gridnum, p, numAtoms, mol, forces, energy, extForce, extVirial);
00146         } else if (grid->get_grid_type() == GridforceGrid::GridforceGridTypeLite) {
00147             GridforceLiteGrid *g = (GridforceLiteGrid *)grid;
00148             do_calc(g, gridnum, p, numAtoms, mol, forces, energy, extForce, extVirial);
00149         }
00150     }
00151     reduction->item(REDUCTION_MISC_ENERGY) += energy;
00152     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
00153     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
00154     reduction->submit();
00155 }
00156 /*                      END OF FUNCTION force                           */

Generated on Thu Sep 21 01:17:10 2017 for NAMD by  doxygen 1.4.7