ComputeEField.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeEField.h"
00009 #include "Node.h"
00010 #include "SimParameters.h"
00011 #include "HomePatch.h"
00012 
00013 
00014 ComputeEField::ComputeEField(ComputeID c, PatchID pid)
00015   : ComputeHomePatch(c,pid)
00016 {
00017 
00018         reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00019 
00020 }
00021 /*                      END OF FUNCTION ComputeEField           */
00022 
00023 
00024 ComputeEField::~ComputeEField()
00025 
00026 {
00027         delete reduction;
00028 }
00029 /*                      END OF FUNCTION ~ComputeEField          */
00030 
00031 
00032 void ComputeEField::doForce(FullAtom* p, Results* r) {
00033 
00034   SimParameters *simParams = Node::Object()->simParameters;
00035   Vector eField = simParams->eField;
00036   // Calculate the angular frequency in 1/fs.
00037   BigReal omega = TWOPI * simParams->eFieldFreq / 1000.;
00038   BigReal phi = PI/180.* simParams->eFieldPhase;
00039   BigReal t = patch->flags.step * simParams->dt;
00040   Vector eField1 = cos(omega * t - phi) * eField;
00041 
00042   const int normalized = simParams->eFieldNormalized;
00043   if ( normalized ) {
00044     Lattice &l = homePatch->lattice;
00045     eField1 = Vector(l.a_r()*eField1, l.b_r()*eField1, l.c_r()*eField1);
00046   }
00047 
00048   Force *forces = r->f[Results::normal];
00049   BigReal energy = 0;
00050   Force extForce = 0.;
00051   Tensor extVirial;
00052 
00053   //  Loop through and check each atom
00054   for (int i=0; i<numAtoms; i++) {
00055     Force force = p[i].charge * eField1; 
00056     forces[i] += force;
00057     Position vpos = homePatch->lattice.reverse_transform(
00058                 p[i].position, p[i].transform );
00059     energy -= force * (vpos - homePatch->lattice.origin());
00060     if ( ! normalized ) {
00061       extForce += force;
00062       extVirial += outer(force,vpos);
00063     }
00064   }
00065 
00066   reduction->item(REDUCTION_MISC_ENERGY) += energy;
00067   if ( ! normalized ) {
00068     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
00069     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
00070   }
00071   reduction->submit();
00072 
00073 }
00074 /*                      END OF FUNCTION force                           */

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