GlobalMasterSMD.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "Node.h"
00009 #include "Molecule.h"
00010 #include "NamdTypes.h"
00011 
00012 #include "GlobalMaster.h"
00013 #include "GlobalMasterSMD.h"
00014 #include "PDB.h"
00015 #include "PDBData.h"
00016 
00017 /* XXX necessary?
00018 #include <iostream.h>
00019 #if !defined(WIN32) || defined(__CYGWIN__)
00020 #include <strstream.h>
00021 #else
00022 #include <strstrea.h>
00023 #endif
00024 #include <string.h>
00025 */
00026 
00027 //#define DEBUGM
00028 #define MIN_DEBUG_LEVEL 1
00029 #include "Debug.h"
00030 
00031 GlobalMasterSMD::GlobalMasterSMD(BigReal spring_constant,
00032                BigReal transverse_spring_constant, BigReal velocity,
00033                const Vector direction, int output_frequency,
00034                int first_timestep, const char *filename,
00035                int numAtoms) {
00036   DebugM(3,"initialize called\n");
00037   moveVel = velocity;
00038   moveDir = direction;
00039   outputFreq = output_frequency;
00040   k = spring_constant;
00041   k2 = transverse_spring_constant;
00042   currentTime = first_timestep;
00043 
00044   parseAtoms(filename,numAtoms);
00045   iout << iINFO << requestedGroups()[0].size() << " SMD ATOMS\n" << endi;
00046   DebugM(1,"done with initialize\n");
00047 }
00048 
00049 void GlobalMasterSMD::parseAtoms(const char *file, int numTotalAtoms) {
00050   DebugM(3,"parseAtoms called\n");
00051   PDB smdpdb(file);
00052   int numatoms = smdpdb.num_atoms();
00053   if (numatoms < 1) 
00054     NAMD_die("No atoms found in SMDFile\n");
00055   if (numatoms != numTotalAtoms)
00056     NAMD_die("The number of atoms in SMDFile must be equal to the total number of atoms in the structure!");
00057 
00058   // add a single group
00059   modifyRequestedGroups().resize(1);
00060 
00061   // add all the required atoms to the group
00062 
00063   // also compute the center of mass of these atoms
00064   cm.x = cm.y = cm.z = 0;
00065   Molecule *mol = Node::Object()->molecule; // to get masses
00066   BigReal imass = 0; 
00067 
00068   for (int i=0; i<numatoms; i++) {
00069 #ifdef MEM_OPT_VERSION
00070     PDBCoreData *atom = smdpdb.atom(i);
00071 #else
00072     PDBAtom *atom = smdpdb.atom(i); // get an atom from the file
00073 #endif    
00074     if (atom->occupancy()) { // if occupancy is not 0, then add it!
00075       // add the atom to the list
00076       modifyRequestedGroups()[0].add(i);
00077 
00078       // compute the center of mass
00079       BigReal mass = mol->atommass(i); 
00080       cm.x += atom->xcoor()*mass;
00081       cm.y += atom->ycoor()*mass;
00082       cm.z += atom->zcoor()*mass;
00083       imass += mass; 
00084     }
00085   }
00086   if (imass == 0) // we didn't find any!
00087     NAMD_die("SMDFile contained no SMD atoms (atoms w/ nonzero occupancy)\n");
00088 
00089   // now compute the center of mass
00090   cm /= imass;
00091   DebugM(1,"done with parseAtoms\n");
00092 }
00093 
00094 GlobalMasterSMD::~GlobalMasterSMD() { }
00095 
00096 void GlobalMasterSMD::calculate() {
00097   DebugM(3,"calculate called\n");
00098 
00099   // make extra-sure we have been set up correctly
00100   if(!requestedGroups().size() == 1)
00101     NAMD_die("Internal error: uninitialized!");
00102 
00103   // set the correct forces
00104   Position curcm = *getGroupPositionBegin(); // get the center of mass
00105   DebugM(1,"Current CM "<<cm<<"\n");
00106   BigReal diff = (curcm - cm)*moveDir;
00107   // second term below is along transverse direction: -diff*moveDir + (diff*moveDir - (curcm-cm)) = -(curcm-cm)
00108   // so if k = k2 and moveVel = 0, we see that f = -k * (curcm - cm), the desired result
00109   Force f = k*(moveVel*currentTime - diff)*moveDir + k2*(diff*moveDir - (curcm - cm));
00110   modifyGroupForces().resize(1);
00111   modifyGroupForces()[0] = f;
00112 
00113   // print some output sometimes
00114   if (currentTime % outputFreq == 0)
00115     output(currentTime, curcm, f);
00116 
00117   // keep track of the number of timesteps
00118   currentTime++;
00119 
00120   DebugM(1,"done with calculate: force: "<<f<<"\n");
00121 }
00122 
00123 void GlobalMasterSMD::output(int t, Position p, Force f) {
00124   DebugM(3,"output called\n");
00125 
00126   if (t % (100*outputFreq) == 0) 
00127     iout << "SMDTITLE: TS   CURRENT_POSITION         FORCE\n" << endi; 
00128   iout << "SMD  " << t << ' ' << p << ' ' << f*PNPERKCALMOL << '\n' << endi;
00129 
00130   DebugM(1,"done with output\n");
00131 }

Generated on Tue Sep 19 01:17:12 2017 for NAMD by  doxygen 1.4.7