/* * Copyright (C) 2005 by David J. Hardy. All rights reserved. * * tinysim.c * * Demonstrates use of mdio, force, and step libraries. * Very simplified, hard-coded config options (in default.h), * no error checking. */ #include #include #include "mdapi/mdtypes.h" #include "mdio/param.h" #include "mdio/topo.h" #include "mdio/pdbcoord.h" #include "force/force.h" #include "step/step.h" #include "debug/debug.h" #include "tinysim/default.h" 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 */ } int main() { /* file I/O */ mdio_Param *param; mdio_Topo *topo; mdio_Pdbcoord *pdb; /* force evaluation */ Force feval; ForceParam fp; ForceEnergy fe; ForceResult fr; /* time stepping */ Step step; StepParam sp; StepSystem sys; /* need some arrays */ MD_Dvec *pos, *vel, *f; int32 natoms, n, i; /* IMPORTANT: clear the memory for helper objects */ memset(&fp, 0, sizeof(ForceParam)); memset(&fe, 0, sizeof(ForceEnergy)); memset(&fr, 0, sizeof(ForceResult)); memset(&sp, 0, sizeof(StepParam)); memset(&sys, 0, sizeof(StepSystem)); /* read force field parameter and topology files */ param = mdio_createParam(); mdio_readParam(param, PARAM_FILE); topo = mdio_createTopo(); mdio_readTopo(topo, PSF_FILE); mdio_indexParamTopo(topo, param); /* have topo arrays index param arrays */ /* set force field parameter data for ForceParam helper object */ 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; /* set topology data for ForceParam helper object */ 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; /* setup other 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 */ /* read initial positions */ pdb = mdio_createPdbcoord(); mdio_readPdbcoord(pdb, PDB_FILE, natoms); /* obtain position array */ pos = mdio_getPdbcoord(pdb, &n); ASSERT(n == natoms); /* must allocate memory for velocity and force arrays */ vel = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec)); f = (MD_Dvec *) calloc(natoms, sizeof(MD_Dvec)); /* setup ForceResult object */ fr.f = f; /* initialize Force class */ force_init(&feval); force_setup(&feval, &fp, &fe, &fr); /* setup StepParam object */ 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.atom = mdio_getAtomTopo(topo, &n); sp.natoms = n; sp.ndegfreedom = 3*natoms - 3; /* subtract 3 because we will be removing center of mass motion */ /* setup StepSystem object */ sys.force = f; sys.vel = vel; sys.pos = pos; /* initialize Step class */ step_init(&step, &sp); /* determine initial velocities */ step_set_random_velocities(&step, &sys, TEMPERATURE); step_remove_com_motion(&step, &sys); /* perform time stepping */ for (i = 0; i < NUMSTEPS; i++) { /* print results after each step */ step_compute(&step, &sys, 1); step_find_reductions(&step, &sys); if (i % LABELSTEPS == 0) { /* print label line */ printf("#%5s %11s %11s %11s %11s\n", "step", "energy", "temp", "ke", "pe"); } 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); } /* 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); return 0; }