NAMD
GlobalMasterSMD.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "Node.h"
9 #include "Molecule.h"
10 #include "NamdTypes.h"
11 
12 #include "GlobalMaster.h"
13 #include "GlobalMasterSMD.h"
14 #include "PDB.h"
15 #include "PDBData.h"
16 
17 /* XXX necessary?
18 #include <iostream.h>
19 #if !defined(WIN32) || defined(__CYGWIN__)
20 #include <strstream.h>
21 #else
22 #include <strstrea.h>
23 #endif
24 #include <string.h>
25 */
26 
27 //#define DEBUGM
28 #define MIN_DEBUG_LEVEL 1
29 #include "Debug.h"
30 
32  BigReal transverse_spring_constant, BigReal velocity,
33  const Vector direction, int output_frequency,
34  int first_timestep, const char *filename,
35  int numAtoms) {
36  DebugM(3,"initialize called\n");
37  moveVel = velocity;
38  moveDir = direction;
39  outputFreq = output_frequency;
40  k = spring_constant;
41  k2 = transverse_spring_constant;
42  currentTime = first_timestep;
43 
44  parseAtoms(filename,numAtoms);
45  iout << iINFO << requestedGroups()[0].size() << " SMD ATOMS\n" << endi;
46  DebugM(1,"done with initialize\n");
47 }
48 
49 void GlobalMasterSMD::parseAtoms(const char *file, int numTotalAtoms) {
50  DebugM(3,"parseAtoms called\n");
51  PDB smdpdb(file);
52  int numatoms = smdpdb.num_atoms();
53  if (numatoms < 1)
54  NAMD_die("No atoms found in SMDFile\n");
55  if (numatoms != numTotalAtoms)
56  NAMD_die("The number of atoms in SMDFile must be equal to the total number of atoms in the structure!");
57 
58  // add a single group
59  modifyRequestedGroups().resize(1);
60 
61  // add all the required atoms to the group
62 
63  // also compute the center of mass of these atoms
64  cm.x = cm.y = cm.z = 0;
65  Molecule *mol = Node::Object()->molecule; // to get masses
66  BigReal imass = 0;
67 
68  for (int i=0; i<numatoms; i++) {
69 #ifdef MEM_OPT_VERSION
70  PDBCoreData *atom = smdpdb.atom(i);
71 #else
72  PDBAtom *atom = smdpdb.atom(i); // get an atom from the file
73 #endif
74  if (atom->occupancy()) { // if occupancy is not 0, then add it!
75  // add the atom to the list
76  modifyRequestedGroups()[0].add(i);
77 
78  // compute the center of mass
79  BigReal mass = mol->atommass(i);
80  cm.x += atom->xcoor()*mass;
81  cm.y += atom->ycoor()*mass;
82  cm.z += atom->zcoor()*mass;
83  imass += mass;
84  }
85  }
86  if (imass == 0) // we didn't find any!
87  NAMD_die("SMDFile contained no SMD atoms (atoms w/ nonzero occupancy)\n");
88 
89  // now compute the center of mass
90  cm /= imass;
91  DebugM(1,"done with parseAtoms\n");
92 }
93 
95 
96 void GlobalMasterSMD::calculate() {
97  DebugM(3,"calculate called\n");
98 
99  // make extra-sure we have been set up correctly
100  if(!requestedGroups().size() == 1)
101  NAMD_die("Internal error: uninitialized!");
102 
103  // set the correct forces
104  Position curcm = *getGroupPositionBegin(); // get the center of mass
105  DebugM(1,"Current CM "<<cm<<"\n");
106  BigReal diff = (curcm - cm)*moveDir;
107  // second term below is along transverse direction: -diff*moveDir + (diff*moveDir - (curcm-cm)) = -(curcm-cm)
108  // so if k = k2 and moveVel = 0, we see that f = -k * (curcm - cm), the desired result
109  Force f = k*(moveVel*currentTime - diff)*moveDir + k2*(diff*moveDir - (curcm - cm));
111  modifyGroupForces()[0] = f;
112 
113  // print some output sometimes
114  if (currentTime % outputFreq == 0)
115  output(currentTime, curcm, f);
116 
117  // keep track of the number of timesteps
118  currentTime++;
119 
120  DebugM(1,"done with calculate: force: "<<f<<"\n");
121 }
122 
123 void GlobalMasterSMD::output(int t, Position p, Force f) {
124  DebugM(3,"output called\n");
125 
126  if (t % (100*outputFreq) == 0)
127  iout << "SMDTITLE: TS CURRENT_POSITION FORCE\n" << endi;
128  iout << "SMD " << t << ' ' << p << ' ' << f*PNPERKCALMOL << '\n' << endi;
129 
130  DebugM(1,"done with output\n");
131 }
static Node * Object()
Definition: Node.h:86
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Definition: PDB.h:35
Definition: Vector.h:64
static int numatoms
Definition: ScriptTcl.C:64
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
BigReal zcoor(void)
Definition: PDBData.C:433
#define iout
Definition: InfoStream.h:87
const ResizeArray< AtomIDList > & requestedGroups()
Definition: GlobalMaster.C:149
GlobalMasterSMD(BigReal spring_constant, BigReal transverse_spring_constant, BigReal velocity, const Vector direction, int output_frequency, int first_timestep, const char *filename, int)
BigReal ycoor(void)
Definition: PDBData.C:429
BigReal x
Definition: Vector.h:66
ResizeArray< AtomIDList > & modifyRequestedGroups()
Definition: GlobalMaster.C:184
void NAMD_die(const char *err_msg)
Definition: common.C:83
ForceList & modifyGroupForces()
Definition: GlobalMaster.C:167
BigReal xcoor(void)
Definition: PDBData.C:425
#define PNPERKCALMOL
Definition: common.h:52
void resize(int i)
Definition: ResizeArray.h:84
BigReal y
Definition: Vector.h:66
Real atommass(int anum) const
Definition: Molecule.h:1038
BigReal occupancy(void)
Definition: PDBData.C:444
infostream & endi(infostream &s)
Definition: InfoStream.C:38
PositionList::const_iterator getGroupPositionBegin()
Definition: GlobalMaster.C:202
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114