MDX Tutorial

previous page / next page

Using the MDIO Library

The MDIO library reads X-Plor force field parameter files and a PSF file, and will then assemble the array representation for the molecular system. There are also routines to read and write PDB coordinate files and NAMD-style binary restart files and routines to write DCD files.

Reading force field parameters:

This is done using the mdio_Param class. 

  mdio_Param *param;

  param = mdio_createParam();
  mdio_readParam(param, PARAM_FILE);

Reading topology:

This is done using the mdio_Topo class. 

  mdio_Topo *topo;

  topo = mdio_createTopo();
  mdio_readTopo(topo, PSF_FILE);

Indexing force field parameters with topology:

The following is considered to be a mdio_Topo class method that also acts on a mdio_Param object. 

  mdio_indexParamTopo(topo, param);

Reading coordinate PDB file:

This is done using the mdio_Pdbcoord class. 

  mdio_Pdbcoord *pdb;

  pdb = mdio_createPdbcoord();
  mdio_readPdbcoord(pdb, PDB_FILE, natoms);  /* natoms is number of atoms */

Using the Force Library

The Force library evaluates a CHARMM-style force field with cutoff nonbonded interactions. You configure its behavior using bit flags to select which potential functions will be contributing to the array of total force. You can also receive the separated force contributions for each atom, if desired. The primary class is Force with private members.

There are three helper classes:

  • ForceParam - used for setup
  • ForceEnergy - provides container for computed potential energies
  • ForceResult - provides container for array of total forces and any desired partial force contributions

Before initializing the helper objects, it is very important to clear their memory contents to all zeros: 

  ForceParam fp;
  ForceEnergy fe;
  ForceResult fr;

  memset(&fp, 0, sizeof(ForceParam));
  memset(&fe, 0, sizeof(ForceEnergy));
  memset(&fr, 0, sizeof(ForceResult));

Use the mdio_Param object to set up the ForceParam force field parameters: 

  fp.atomprm = mdio_getAtomParam(param, &n);
  fp.atomprm_len = n;
  fp.bondprm = mdio_getBondParam(param, &n);
  fp.bondprm_len = n;
  fp.angleprm = mdio_getAngleParam(param, &n);
  fp.angleprm_len = n;
  fp.dihedprm = mdio_getDihedParam(param, &n);
  fp.dihedprm_len = n;
  fp.imprprm = mdio_getImprParam(param, &n);
  fp.imprprm_len = n;
  fp.nbfixprm = mdio_getNbfixParam(param, &n);
  fp.nbfixprm_len = n;

Use the mdio_Topo object to set up the ForceParam topology: 

  fp.atom = mdio_getAtomTopo(topo, &natoms);
  fp.atom_len = natoms;  /* number of atoms in system */
  fp.bond = mdio_getBondTopo(topo, &n);
  fp.bond_len = n;
  fp.angle = mdio_getAngleTopo(topo, &n);
  fp.angle_len = n;
  fp.dihed = mdio_getDihedTopo(topo, &n);
  fp.dihed_len = n;
  fp.impr = mdio_getImprTopo(topo, &n);
  fp.impr_len = n;
  fp.excl = mdio_getExclTopo(topo, &n);
  fp.excl_len = n;

Set up the other needed ForceParam fields: 

  fp.cutoff = CUTOFF;
  fp.switchdist = SWITCHDIST;
  fp.elec_const = MD_COULOMB;  /* MD_COULOMB is from mdtypes.h */
  fp.dielectric = DIELECTRIC;
  fp.scaling14 = SCALING14;
  fp.flags = FORCE_MD_VACUUM;  /* FORCE_MD_VACUUM is collection of flags for
                                  smooth cutoff MD, scaled 1-4 exclusions,
                                  no boundary restraints */

Provide to ForceResult an array buffer space for the total forces: 

  MD_Dvec *f;

  f = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec));
  fr.f = f;

Now that we have the helper objects properly initialized, we can call the Force constructor and setup functions: 

  Force feval;

  force_init(&feval);
  force_setup(&feval, &fp, &fe, &fr);

Using the Step Library

The Step library performs time stepping with leapfrog (velocity Verlet) integration. It also can initialize velocities and compute kinetic energy and temperature reductions. The primary class is Step with private members.

There are two helper classes:

  • StepParam - used for setup
  • StepSystem - provides arrays for trajectories, returns kinetic energy and temperature values

As before, it is very important to reset the memory contents of the helper objects: 

  StepParam sp;
  StepSystem sys;

  memset(&sp, 0, sizeof(StepParam));
  memset(&sys, 0, sizeof(StepSystem));

Use the mdio_Topo object to set up the StepParam topology: 

  sp.atom = mdio_getAtomTopo(topo, &n);
  sp.natoms = n;

Use the mdio_Pdbcoord object to set up the StepSystem force field parameters: 

  MD_Dvec *pos;

  pos = mdio_getPdbcoord(pdb, &n);
  sys.pos = pos;

Provide to StepSystem an array buffer space for the atom velocities: 

  MD_Dvec *vel;

  vel = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec));
  sys.vel = vel;

For this example, the Force class is our only resource for force evaluation, so it is OK to also assign to StepSystem the array buffer space allocated previously for the ForceResult total forces: 

  sys.force = f;

Set up the remaining StepParam fields: 

  sp.random_seed = SEED;
  sp.timestep = TIMESTEP;
  sp.force_compute = my_force_compute;  /* our own force computation routine */
  sp.force_object = (void *) (&feval);  /* and data that must be passed to it */
  sp.ndegfreedom = 3*natoms - 3;        /* subtract 3 because we will be
                                           removing center of mass motion */

This is the tricky part! The my_force_compute from above is a function that is called by the Step method step_compute() in order to receive the system forces. Here is the definition of my_force_compute: 

int my_force_compute(void *vfeval, MD_Dvec *f, MD_Dvec *pos)
{
  Force *feval = (Force *) vfeval;

  /* Force class already has pointer to the f array */
  return force_compute(feval, pos, NULL);  /* set "wrap" array to NULL */
}
We recover the Force object by typecasting the address assigned to sp.force_object and use this to invoke the force_compute() method. Since we previously assigned sys.force to the same array buffer space as fr.f, we do not need to explicitly deal with the force array f expected by step_compute() (it is automatically populated).

Note that if we were doing more than this for the force evaluation, e.g. computing long-range electrostatic forces using the Mgrid multilevel summation solver, we would use my_force_compute() to manage calls to other force evaluation libraries and explicitly sum the results into the f array.

Now we are ready to initialize the Step object and determine the initial velocities: 

  Step step;

  /* initialize Step class */
  step_init(&step, &sp);

  /* determine initial velocities */
  step_set_random_velocities(&step, &sys, TEMPERATURE);
  step_remove_com_motion(&step, &sys);

The time stepping and the reporting of results after each step is done by the following loop: 

  /* perform time stepping */
  for (i = 0;  i < NUMSTEPS;  i++) {

    /* print results after each step */
    step_compute(&step, &sys, 1);
    step_find_reductions(&step, &sys);
    printf(" %5d %11.6f %11.6f %11.6f %11.6f\n", i+1,
        fe.pe + sys.kinetic_energy,
        sys.temperature, sys.kinetic_energy, fe.pe);
  }

Once the simulation is finished, our explicit memory allocations along with that done implicitly through the objects can all be freed: 

  /* free memory our memory allocations */
  free(vel);
  free(f);

  /* finished with objects */
  step_done(&step);
  force_done(&feval);
  mdio_destroyPdbcoord(pdb);
  mdio_destroyTopo(topo);
  mdio_destroyParam(param);

previous page / next page