fix_cmap.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Jul 10, 00:03

fix_cmap.cpp
View Options

	/* ----------------------------------------------------------------------
	LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
	http://lammps.sandia.gov, Sandia National Laboratories
	Steve Plimpton, sjplimp@sandia.gov

	Copyright (2003) Sandia Corporation. Under the terms of Contract
	DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
	certain rights in this software. This software is distributed under
	the GNU General Public License.

	See the README file in the top-level LAMMPS directory.
	------------------------------------------------------------------------- */

	/* ----------------------------------------------------------------------
	Contributing authors:
	Xiaohu Hu, CMB/ORNL (hux2@ornl.gov)
	David Hyde-Volpe, Tigran Abramyan, and Robert A. Latour (Clemson University)
	Chris Lorenz (Kings College-London)

	Implementation of the CHARMM CMAP; adds an extra energy term for the
	peptide backbone dihedrals. The tools/ch2lmp/charmm2lammps.pl
	conversion script, which generates an extra section in the LAMMPS data
	file, is needed in order to generate the info used by this fix style.

	References:
	- MacKerell et al., J. Am. Chem. Soc. 126(2004):698-699.
	- MacKerell et al., J. Comput. Chem. 25(2004):1400-1415.
	------------------------------------------------------------------------- */

	#include <mpi.h>
	#include <math.h>
	#include <stdlib.h>
	#include <string.h>
	#include <stdio.h>
	#include "fix_cmap.h"
	#include "atom.h"
	#include "atom_vec.h"
	#include "update.h"
	#include "respa.h"
	#include "modify.h"
	#include "domain.h"
	#include "force.h"
	#include "group.h"
	#include "comm.h"
	#include "math_const.h"
	#include "memory.h"
	#include "error.h"

	using namespace LAMMPS_NS;
	using namespace FixConst;
	using namespace MathConst;

	#define MAXLINE 256
	#define LISTDELTA 10000
	#define LB_FACTOR 1.5

	#define CMAPMAX 6 // max # of CMAP terms stored by one atom
	#define CMAPDIM 24 // grid map dimension is 24 x 24
	#define CMAPXMIN -360.0
	#define CMAPXMIN2 -180.0
	#define CMAPDX 15.0 // 360/CMAPDIM

	/* ---------------------------------------------------------------------- */

	FixCMAP::FixCMAP(LAMMPS lmp, int narg, char *arg) :
	Fix(lmp, narg, arg),
	crosstermlist(NULL), num_crossterm(NULL), crossterm_type(NULL),
	crossterm_atom1(NULL), crossterm_atom2(NULL), crossterm_atom3(NULL),
	crossterm_atom4(NULL), crossterm_atom5(NULL),
	g_axis(NULL), cmapgrid(NULL), d1cmapgrid(NULL), d2cmapgrid(NULL),
	d12cmapgrid(NULL)
	{
	if (narg != 4) error->all(FLERR,"Illegal fix cmap command");

	restart_global = 1;
	restart_peratom = 1;
	peatom_flag = 1;
	virial_flag = 1;
	peratom_freq = 1;
	scalar_flag = 1;
	global_freq = 1;
	extscalar = 1;
	extvector = 1;
	wd_header = 1;
	wd_section = 1;

	MPI_Comm_rank(world,&me);
	MPI_Comm_size(world,&nprocs);

	// allocate memory for CMAP data

	memory->create(g_axis,CMAPDIM,"cmap:g_axis");
	memory->create(cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:grid");
	memory->create(d1cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d1grid");
	memory->create(d2cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d2grid");
	memory->create(d12cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d12grid");

	// read and setup CMAP data

	read_grid_map(arg[3]);

	// perform initial allocation of atom-based arrays
	// register with Atom class

	num_crossterm = NULL;
	crossterm_type = NULL;
	crossterm_atom1 = NULL;
	crossterm_atom2 = NULL;
	crossterm_atom3 = NULL;
	crossterm_atom4 = NULL;
	crossterm_atom5 = NULL;

	nmax_previous = 0;
	grow_arrays(atom->nmax);
	atom->add_callback(0);
	atom->add_callback(1);

	// local list of crossterms

	ncmap = 0;
	maxcrossterm = 0;
	crosstermlist = NULL;
	}

	/* --------------------------------------------------------------------- */

	FixCMAP::~FixCMAP()
	{
	// unregister callbacks to this fix from Atom class

	atom->delete_callback(id,0);
	atom->delete_callback(id,1);

	memory->destroy(g_axis);
	memory->destroy(cmapgrid);
	memory->destroy(d1cmapgrid);
	memory->destroy(d2cmapgrid);
	memory->destroy(d12cmapgrid);

	memory->destroy(crosstermlist);

	memory->destroy(num_crossterm);
	memory->destroy(crossterm_type);
	memory->destroy(crossterm_atom1);
	memory->destroy(crossterm_atom2);
	memory->destroy(crossterm_atom3);
	memory->destroy(crossterm_atom4);
	memory->destroy(crossterm_atom5);
	}

	/* ---------------------------------------------------------------------- */

	int FixCMAP::setmask()
	{
	int mask = 0;
	mask \|= PRE_NEIGHBOR;
	mask \|= PRE_REVERSE;
	mask \|= POST_FORCE;
	mask \|= THERMO_ENERGY;
	mask \|= POST_FORCE_RESPA;
	mask \|= MIN_POST_FORCE;
	return mask;
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::init()
	{
	int i;
	double angle;

	i = 0;
	angle = -180.0;
	while (angle < 180.0) {
	g_axis[i] = angle;
	angle += CMAPDX;
	i++;
	}

	// pre-compute the derivatives of the maps

	for (i = 0; i < 6; i++)
	set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]);

	// define newton_bond here in case restart file was read (not data file)

	newton_bond = force->newton_bond;
	}

	/* --------------------------------------------------------------------- */

	void FixCMAP::setup(int vflag)
	{
	pre_neighbor();

	if (strstr(update->integrate_style,"verlet"))
	post_force(vflag);
	else {
	((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
	post_force_respa(vflag,nlevels_respa-1,0);
	((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
	}
	}

	/* --------------------------------------------------------------------- */

	void FixCMAP::setup_pre_neighbor()
	{
	pre_neighbor();
	}

	/* --------------------------------------------------------------------- */

	void FixCMAP::setup_pre_reverse(int eflag, int vflag)
	{
	pre_reverse(eflag,vflag);
	}

	/* --------------------------------------------------------------------- */

	void FixCMAP::min_setup(int vflag)
	{
	pre_neighbor();
	post_force(vflag);
	}

	/* ----------------------------------------------------------------------
	store local neighbor list as if newton_bond = OFF, even if actually ON
	------------------------------------------------------------------------- */

	void FixCMAP::pre_neighbor()
	{
	int i,m,atom1,atom2,atom3,atom4,atom5;

	// guesstimate initial length of local crossterm list
	// if ncmap was not set (due to read_restart, no read_data),
	// then list will grow by LISTDELTA chunks

	if (maxcrossterm == 0) {
	if (nprocs == 1) maxcrossterm = ncmap;
	else maxcrossterm = static_cast<int> (LB_FACTOR*ncmap/nprocs);
	memory->create(crosstermlist,maxcrossterm,6,"cmap:crosstermlist");
	}

	int nlocal = atom->nlocal;

	ncrosstermlist = 0;

	for (i = 0; i < nlocal; i++) {
	for (m = 0; m < num_crossterm[i]; m++) {
	atom1 = atom->map(crossterm_atom1[i][m]);
	atom2 = atom->map(crossterm_atom2[i][m]);
	atom3 = atom->map(crossterm_atom3[i][m]);
	atom4 = atom->map(crossterm_atom4[i][m]);
	atom5 = atom->map(crossterm_atom5[i][m]);

	if (atom1 == -1 \|\| atom2 == -1 \|\| atom3 == -1 \|\|
	atom4 == -1 \|\| atom5 == -1) {
	char str[128];
	sprintf(str,"CMAP atoms "
	TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
	TAGINT_FORMAT " " TAGINT_FORMAT
	" missing on proc %d at step " BIGINT_FORMAT,
	crossterm_atom1[i][m],crossterm_atom2[i][m],
	crossterm_atom3[i][m],crossterm_atom4[i][m],
	crossterm_atom5[i][m],me,update->ntimestep);
	error->one(FLERR,str);
	}
	atom1 = domain->closest_image(i,atom1);
	atom2 = domain->closest_image(i,atom2);
	atom3 = domain->closest_image(i,atom3);
	atom4 = domain->closest_image(i,atom4);
	atom5 = domain->closest_image(i,atom5);

	if (i <= atom1 && i <= atom2 && i <= atom3 &&
	i <= atom4 && i <= atom5) {
	if (ncrosstermlist == maxcrossterm) {
	maxcrossterm += LISTDELTA;
	memory->grow(crosstermlist,maxcrossterm,6,"cmap:crosstermlist");
	}
	crosstermlist[ncrosstermlist][0] = atom1;
	crosstermlist[ncrosstermlist][1] = atom2;
	crosstermlist[ncrosstermlist][2] = atom3;
	crosstermlist[ncrosstermlist][3] = atom4;
	crosstermlist[ncrosstermlist][4] = atom5;
	crosstermlist[ncrosstermlist][5] = crossterm_type[i][m];
	ncrosstermlist++;
	}
	}
	}
	}

	/* ----------------------------------------------------------------------
	store eflag, so can use it in post_force to tally per-atom energies
	------------------------------------------------------------------------- */

	void FixCMAP::pre_reverse(int eflag, int vflag)
	{
	eflag_caller = eflag;
	}

	/* ----------------------------------------------------------------------
	compute CMAP terms as if newton_bond = OFF, even if actually ON
	------------------------------------------------------------------------- */

	void FixCMAP::post_force(int vflag)
	{
	int n,i1,i2,i3,i4,i5,type,nlist;
	int li1, li2, mli1,mli2,mli11,mli21,t1,li3,li4,mli3,mli4,mli31,mli41;
	int list[5];
	// vectors needed to calculate the cross-term dihedral angles
	double vb21x,vb21y,vb21z,vb32x,vb32y,vb32z,vb34x,vb34y,vb34z;
	double vb23x,vb23y,vb23z;
	double vb43x,vb43y,vb43z,vb45x,vb45y,vb45z,a1x,a1y,a1z,b1x,b1y,b1z;
	double a2x,a2y,a2z,b2x,b2y,b2z,r32,a1sq,b1sq,a2sq,b2sq,dpr21r32,dpr34r32;
	double dpr32r43,dpr45r43,r43,vb12x,vb12y,vb12z,vb54x,vb54y,vb54z;
	// cross-term dihedral angles
	double phi,psi,phi1,psi1;
	double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6];
	double gs[4],d1gs[4],d2gs[4],d12gs[4];
	double engfraction;
	// vectors needed for the gradient/force calculation
	double dphidr1x,dphidr1y,dphidr1z,dphidr2x,dphidr2y,dphidr2z;
	double dphidr3x,dphidr3y,dphidr3z,dphidr4x,dphidr4y,dphidr4z;
	double dpsidr1x,dpsidr1y,dpsidr1z,dpsidr2x,dpsidr2y,dpsidr2z;
	double dpsidr3x,dpsidr3y,dpsidr3z,dpsidr4x,dpsidr4y,dpsidr4z;

	// Definition of cross-term dihedrals

	// phi dihedral
	// \|--------------------\|
	// a1-----a2-----a3-----a4-----a5 cross-term atoms
	// C N CA C N cross-term atom types
	// \|--------------------\|
	// psi dihedral

	double **x = atom->x;
	double **f = atom->f;
	int nlocal = atom->nlocal;

	ecmap = 0.0;
	int eflag = eflag_caller;
	if (eflag \|\| vflag) ev_setup(eflag,vflag);
	else evflag = 0;

	for (n = 0; n < ncrosstermlist; n++) {
	i1 = crosstermlist[n][0];
	i2 = crosstermlist[n][1];
	i3 = crosstermlist[n][2];
	i4 = crosstermlist[n][3];
	i5 = crosstermlist[n][4];

	type = crosstermlist[n][5];
	if (type == 0) continue;

	// calculate bond vectors for both dihedrals

	// phi
	// vb21 = r2 - r1

	vb21x = x[i2][0] - x[i1][0];
	vb21y = x[i2][1] - x[i1][1];
	vb21z = x[i2][2] - x[i1][2];
	vb12x = -1.0*vb21x;
	vb12y = -1.0*vb21y;
	vb12z = -1.0*vb21z;
	vb32x = x[i3][0] - x[i2][0];
	vb32y = x[i3][1] - x[i2][1];
	vb32z = x[i3][2] - x[i2][2];
	vb23x = -1.0*vb32x;
	vb23y = -1.0*vb32y;
	vb23z = -1.0*vb32z;

	vb34x = x[i3][0] - x[i4][0];
	vb34y = x[i3][1] - x[i4][1];
	vb34z = x[i3][2] - x[i4][2];

	// psi
	// bond vectors same as for phi: vb32

	vb43x = -1.0*vb34x;
	vb43y = -1.0*vb34y;
	vb43z = -1.0*vb34z;

	vb45x = x[i4][0] - x[i5][0];
	vb45y = x[i4][1] - x[i5][1];
	vb45z = x[i4][2] - x[i5][2];
	vb54x = -1.0*vb45x;
	vb54y = -1.0*vb45y;
	vb54z = -1.0*vb45z;

	// calculate normal vectors for planes that define the dihedral angles

	a1x = vb12yvb23z - vb12zvb23y;
	a1y = vb12zvb23x - vb12xvb23z;
	a1z = vb12xvb23y - vb12yvb23x;

	b1x = vb43yvb23z - vb43zvb23y;
	b1y = vb43zvb23x - vb43xvb23z;
	b1z = vb43xvb23y - vb43yvb23x;

	a2x = vb23yvb34z - vb23zvb34y;
	a2y = vb23zvb34x - vb23xvb34z;
	a2z = vb23xvb34y - vb23yvb34x;

	b2x = vb45yvb43z - vb45zvb43y;
	b2y = vb45zvb43x - vb45xvb43z;
	b2z = vb45xvb43y - vb45yvb43x;

	// calculate terms used later in calculations

	r32 = sqrt(vb32xvb32x + vb32yvb32y + vb32z*vb32z);
	a1sq = a1xa1x + a1ya1y + a1z*a1z;
	b1sq = b1xb1x + b1yb1y + b1z*b1z;

	r43 = sqrt(vb43xvb43x + vb43yvb43y + vb43z*vb43z);
	a2sq = a2xa2x + a2ya2y + a2z*a2z;
	b2sq = b2xb2x + b2yb2y + b2z*b2z;
	//if (a1sq<0.0001 \|\| b1sq<0.0001 \|\| a2sq<0.0001 \|\| b2sq<0.0001)
	// printf("a1sq b1sq a2sq b2sq: %f %f %f %f \n",a1sq,b1sq,a2sq,b2sq);
	if (a1sq<0.0001 \|\| b1sq<0.0001 \|\| a2sq<0.0001 \|\| b2sq<0.0001) continue;
	dpr21r32 = vb21xvb32x + vb21yvb32y + vb21z*vb32z;
	dpr34r32 = vb34xvb32x + vb34yvb32y + vb34z*vb32z;
	dpr32r43 = vb32xvb43x + vb32yvb43y + vb32z*vb43z;
	dpr45r43 = vb45xvb43x + vb45yvb43y + vb45z*vb43z;

	// calculate the backbone dihedral angles as VMD and GROMACS

	phi = dihedral_angle_atan2(vb21x,vb21y,vb21z,a1x,a1y,a1z,b1x,b1y,b1z,r32);
	psi = dihedral_angle_atan2(vb32x,vb32y,vb32z,a2x,a2y,a2z,b2x,b2y,b2z,r43);

	if (phi == 180.0) phi= -180.0;
	if (psi == 180.0) psi= -180.0;

	phi1 = phi;
	if (phi1 < 0.0) phi1 += 360.0;
	psi1 = psi;
	if (psi1 < 0.0) psi1 += 360.0;

	// find the neighbor grid point index

	li1 = int(((phi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));
	li2 = int(((psi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));

	li3 = int((phi-CMAPXMIN2)/CMAPDX);
	li4 = int((psi-CMAPXMIN2)/CMAPDX);
	mli3 = li3 % CMAPDIM;
	mli4 = li4 % CMAPDIM;
	mli31 = (li3+1) % CMAPDIM;
	mli41 = (li4+1) %CMAPDIM;
	mli1 = li1 % CMAPDIM;
	mli2 = li2 % CMAPDIM;
	mli11 = (li1+1) % CMAPDIM;
	mli21 = (li2+1) %CMAPDIM;
	t1 = type-1;
	if (t1 < 0 \|\| t1 > 5) error->all(FLERR,"Invalid CMAP crossterm_type");

	// determine the values and derivatives for the grid square points

	gs[0] = cmapgrid[t1][mli3][mli4];
	gs[1] = cmapgrid[t1][mli31][mli4];
	gs[2] = cmapgrid[t1][mli31][mli41];
	gs[3] = cmapgrid[t1][mli3][mli41];
	d1gs[0] = d1cmapgrid[t1][mli1][mli2];
	d1gs[1] = d1cmapgrid[t1][mli11][mli2];
	d1gs[2] = d1cmapgrid[t1][mli11][mli21];
	d1gs[3] = d1cmapgrid[t1][mli1][mli21];

	d2gs[0] = d2cmapgrid[t1][mli1][mli2];
	d2gs[1] = d2cmapgrid[t1][mli11][mli2];
	d2gs[2] = d2cmapgrid[t1][mli11][mli21];
	d2gs[3] = d2cmapgrid[t1][mli1][mli21];

	d12gs[0] = d12cmapgrid[t1][mli1][mli2];
	d12gs[1] = d12cmapgrid[t1][mli11][mli2];
	d12gs[2] = d12cmapgrid[t1][mli11][mli21];
	d12gs[3] = d12cmapgrid[t1][mli1][mli21];

	// calculate the cmap energy and the gradient (dE/dphi,dE/dpsi)

	bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs);

	// sum up cmap energy contributions

	engfraction = 0.2 * E;
	if (i1 < nlocal) ecmap += engfraction;
	if (i2 < nlocal) ecmap += engfraction;
	if (i3 < nlocal) ecmap += engfraction;
	if (i4 < nlocal) ecmap += engfraction;
	if (i5 < nlocal) ecmap += engfraction;

	// calculate the derivatives dphi/dr_i

	dphidr1x = 1.0r32/a1sqa1x;
	dphidr1y = 1.0r32/a1sqa1y;
	dphidr1z = 1.0r32/a1sqa1z;

	dphidr2x = -1.0r32/a1sqa1x - dpr21r32/a1sq/r32*a1x +
	dpr34r32/b1sq/r32*b1x;
	dphidr2y = -1.0r32/a1sqa1y - dpr21r32/a1sq/r32*a1y +
	dpr34r32/b1sq/r32*b1y;
	dphidr2z = -1.0r32/a1sqa1z - dpr21r32/a1sq/r32*a1z +
	dpr34r32/b1sq/r32*b1z;

	dphidr3x = dpr34r32/b1sq/r32b1x - dpr21r32/a1sq/r32a1x - r32/b1sq*b1x;
	dphidr3y = dpr34r32/b1sq/r32b1y - dpr21r32/a1sq/r32a1y - r32/b1sq*b1y;
	dphidr3z = dpr34r32/b1sq/r32b1z - dpr21r32/a1sq/r32a1z - r32/b1sq*b1z;

	dphidr4x = r32/b1sq*b1x;
	dphidr4y = r32/b1sq*b1y;
	dphidr4z = r32/b1sq*b1z;

	// calculate the derivatives dpsi/dr_i

	dpsidr1x = 1.0r43/a2sqa2x;
	dpsidr1y = 1.0r43/a2sqa2y;
	dpsidr1z = 1.0r43/a2sqa2z;

	dpsidr2x = r43/a2sqa2x + dpr32r43/a2sq/r43a2x - dpr45r43/b2sq/r43*b2x;
	dpsidr2y = r43/a2sqa2y + dpr32r43/a2sq/r43a2y - dpr45r43/b2sq/r43*b2y;
	dpsidr2z = r43/a2sqa2z + dpr32r43/a2sq/r43a2z - dpr45r43/b2sq/r43*b2z;

	dpsidr3x = dpr45r43/b2sq/r43b2x - dpr32r43/a2sq/r43a2x - r43/b2sq*b2x;
	dpsidr3y = dpr45r43/b2sq/r43b2y - dpr32r43/a2sq/r43a2y - r43/b2sq*b2y;
	dpsidr3z = dpr45r43/b2sq/r43b2z - dpr32r43/a2sq/r43a2z - r43/b2sq*b2z;

	dpsidr4x = r43/b2sq*b2x;
	dpsidr4y = r43/b2sq*b2y;
	dpsidr4z = r43/b2sq*b2z;

	// calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr)

	f1[0] = dEdPhi*dphidr1x;
	f1[1] = dEdPhi*dphidr1y;
	f1[2] = dEdPhi*dphidr1z;
	f2[0] = dEdPhidphidr2x + dEdPsidpsidr1x;
	f2[1] = dEdPhidphidr2y + dEdPsidpsidr1y;
	f2[2] = dEdPhidphidr2z + dEdPsidpsidr1z;
	f3[0] = -dEdPhidphidr3x - dEdPsidpsidr2x;
	f3[1] = -dEdPhidphidr3y - dEdPsidpsidr2y;
	f3[2] = -dEdPhidphidr3z - dEdPsidpsidr2z;
	f4[0] = -dEdPhidphidr4x - dEdPsidpsidr3x;
	f4[1] = -dEdPhidphidr4y - dEdPsidpsidr3y;
	f4[2] = -dEdPhidphidr4z - dEdPsidpsidr3z;
	f5[0] = -dEdPsi*dpsidr4x;
	f5[1] = -dEdPsi*dpsidr4y;
	f5[2] = -dEdPsi*dpsidr4z;

	// apply force to each of the 5 atoms

	if (i1 < nlocal) {
	f[i1][0] += f1[0];
	f[i1][1] += f1[1];
	f[i1][2] += f1[2];
	}
	if (i2 < nlocal) {
	f[i2][0] += f2[0];
	f[i2][1] += f2[1];
	f[i2][2] += f2[2];
	}
	if (i3 < nlocal) {
	f[i3][0] += f3[0];
	f[i3][1] += f3[1];
	f[i3][2] += f3[2];
	}
	if (i4 < nlocal) {
	f[i4][0] += f4[0];
	f[i4][1] += f4[1];
	f[i4][2] += f4[2];
	}
	if (i5 < nlocal) {
	f[i5][0] += f5[0];
	f[i5][1] += f5[1];
	f[i5][2] += f5[2];
	}

	// tally energy and/or virial

	if (evflag) {
	nlist = 0;
	if (i1 < nlocal) list[nlist++] = i1;
	if (i2 < nlocal) list[nlist++] = i2;
	if (i3 < nlocal) list[nlist++] = i3;
	if (i4 < nlocal) list[nlist++] = i4;
	if (i5 < nlocal) list[nlist++] = i5;
	vcmap[0] = (vb12xf1[0])+(vb32xf3[0])+((vb43x+vb32x)*f4[0])+
	((vb54x+vb43x+vb32x)*f5[0]);
	vcmap[1] = (vb12yf1[1])+(vb32yf3[1])+((vb43y+vb32y)*f4[1])+
	((vb54y+vb43y+vb32y)*f5[1]);
	vcmap[2] = (vb12zf1[2])+(vb32zf3[2])+((vb43z+vb32z)*f4[2])+
	((vb54z+vb43z+vb32z)*f5[2]);
	vcmap[3] = (vb12xf1[1])+(vb32xf3[1])+((vb43x+vb32x)*f4[1])+
	((vb54x+vb43x+vb32x)*f5[1]);
	vcmap[4] = (vb12xf1[2])+(vb32xf3[2])+((vb43x+vb32x)*f4[2])+
	((vb54x+vb43x+vb32x)*f5[2]);
	vcmap[5] = (vb12yf1[2])+(vb32yf3[2])+((vb43y+vb32y)*f4[2])+
	((vb54y+vb43y+vb32y)*f5[2]);
	ev_tally(nlist,list,5.0,E,vcmap);
	//ev_tally(5,list,nlocal,newton_bond,E,vcmap);
	}
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::post_force_respa(int vflag, int ilevel, int iloop)
	{
	if (ilevel == nlevels_respa-1) post_force(vflag);
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::min_post_force(int vflag)
	{
	post_force(vflag);
	}

	/* ----------------------------------------------------------------------
	energy of CMAP term
	------------------------------------------------------------------------- */

	double FixCMAP::compute_scalar()
	{
	double all;
	MPI_Allreduce(&ecmap,&all,1,MPI_DOUBLE,MPI_SUM,world);
	return all;
	}

	// ----------------------------------------------------------------------
	// ----------------------------------------------------------------------
	// methods to read CMAP potential file, perform interpolation
	// ----------------------------------------------------------------------
	// ----------------------------------------------------------------------

	void FixCMAP::read_grid_map(char *cmapfile)
	{
	char linebuf[MAXLINE];
	char chunk,line;
	int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter;

	FILE *fp = NULL;
	if (comm->me == 0) {
	fp = force->open_potential(cmapfile);
	if (fp == NULL) {
	char str[128];
	sprintf(str,"Cannot open fix cmap file %s",cmapfile);
	error->one(FLERR,str);
	}
	}

	for (int ix1 = 0; ix1 < 6; ix1++)
	for (int ix2 = 0; ix2 < CMAPDIM; ix2++)
	for (int ix3 = 0; ix3 < CMAPDIM; ix3++)
	cmapgrid[ix1][ix2][ix3] = 0.0;

	counter = 0;
	i1 = i2 = i3 = i4 = i5 = i6 = 0;
	j1 = j2 = j3 = j4 = j5 = j6 = 0;

	int done = 0;

	while (!done) {
	// only read on rank 0 and broadcast to all other ranks
	if (comm->me == 0)
	done = (fgets(linebuf,MAXLINE,fp) == NULL);

	MPI_Bcast(&done,1,MPI_INT,0,world);
	if (done) continue;

	MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world);

	// remove leading whitespace
	line = linebuf;
	while (line && (line == ' ' \|\| line == '\t' \|\| *line == '\r')) ++line;

	// skip if empty line or comment
	if (!line \|\| line =='\n' \|\| line == '\0' \|\| *line == '#') continue;

	// read in the cmap grid point values
	// NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors
	// will occur if content of the file "cmap.data" is altered
	//
	// Reading order of the maps:
	// 1. Alanine map
	// 2. Alanine before proline map
	// 3. Proline map
	// 4. Two adjacent prolines map
	// 5. Glycine map
	// 6. Glycine before proline map

	chunk = strtok(line, " \r\n");
	while (chunk != NULL) {

	// alanine map

	if (counter < CMAPDIM*CMAPDIM) {
	cmapgrid[0][i1][j1] = atof(chunk);
	chunk = strtok(NULL, " \r\n");
	j1++;
	if (j1 == CMAPDIM) {
	j1 = 0;
	i1++;
	}
	counter++;
	}

	// alanine-proline map

	else if (counter >= CMAPDIM*CMAPDIM &&
	counter < 2CMAPDIMCMAPDIM) {
	cmapgrid[1][i2][j2]= atof(chunk);
	chunk = strtok(NULL, " \r\n");
	j2++;
	if (j2 == CMAPDIM) {
	j2 = 0;
	i2++;
	}
	counter++;
	}

	// proline map

	else if (counter >= 2CMAPDIMCMAPDIM &&
	counter < 3CMAPDIMCMAPDIM) {
	cmapgrid[2][i3][j3] = atof(chunk);
	chunk = strtok(NULL, " \r\n");
	j3++;
	if (j3 == CMAPDIM) {
	j3 = 0;
	i3++;
	}
	counter++;
	}

	// 2 adjacent prolines map

	else if (counter >= 3CMAPDIMCMAPDIM &&
	counter < 4CMAPDIMCMAPDIM) {
	cmapgrid[3][i4][j4] = atof(chunk);
	chunk = strtok(NULL, " \r\n");
	j4++;
	if (j4 == CMAPDIM) {
	j4 = 0;
	i4++;
	}
	counter++;
	}

	// glycine map

	else if (counter >= 4CMAPDIMCMAPDIM &&
	counter < 5CMAPDIMCMAPDIM) {
	cmapgrid[4][i5][j5] = atof(chunk);
	chunk = strtok(NULL, " \r\n");
	j5++;
	if (j5 == CMAPDIM) {
	j5 = 0;
	i5++;
	}
	counter++;
	}

	// glycine-proline map

	else if (counter >= 5CMAPDIMCMAPDIM &&
	counter < 6CMAPDIMCMAPDIM) {
	cmapgrid[5][i6][j6] = atof(chunk);
	chunk = strtok(NULL, " \r\n");
	j6++;
	if (j6 == CMAPDIM) {
	j6 = 0;
	i6++;
	}
	counter++;
	}

	else break;
	}
	}

	if (comm->me == 0) fclose(fp);
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::spline(double y, double ddy, int n)
	{
	// create the 2nd dervatives of a taublated function y_i(x_i)
	// at the tabulated points

	int i, j;
	double p, *u;

	memory->create(u,n-1,"cmap:u");

	ddy[0] = u[0] = 0.0;

	for (i = 1; i <= n-2; i++) {
	p = 1.0/(ddy[i-1]+4.0);
	ddy[i] = -p;
	u[i] = ((((6.0y[i+1])-(12.0y[i])+(6.0y[i-1]))/(CMAPDXCMAPDX))-u[i-1])*p;
	}

	ddy[n-1] = 0.0;

	for (j = n-2; j >= 0; j--)
	ddy[j] = ddy[j]*ddy[j+1] + u[j];

	memory->destroy(u);
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::spl_interpolate(double x, double y, double ddy, double &yo,
	double &dyo)
	{
	// perform a 1D cubic spline interpolation

	int ix;
	double a,b,a1,b1,a2,b2;

	ix = int((x-CMAPXMIN)/CMAPDX-(1./2.));

	a = (CMAPXMIN+(ix1.0)CMAPDX-x)/CMAPDX;
	b = (x-CMAPXMIN-(((ix-1)1.0)CMAPDX))/CMAPDX;

	a1 = aaa-a;
	b1 = bbb-b;

	a2 = 3.0aa-1.0;
	b2 = 3.0bb-1.0;
	yo = ay[ix]+by[ix+1]+(a1ddy[ix]+b1ddy[ix+1])(CMAPDXCMAPDX)/6.0;
	dyo = (y[ix+1]-y[ix])/CMAPDX-a2/6.0CMAPDXddy[ix]+b2/6.0CMAPDXddy[ix+1];
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::set_map_derivatives(double map, double d1yo, double **d2yo,
	double **d12yo)
	{
	// precompute the gradient and cross-derivatives of the map grid points.
	// use the bicubic spline to calculate the derivatives

	int i, j, k, ii, jj, xm, p;
	double phi, psi, d1y, d2y, d12y, tyyk,tdyk;
	double tmp_y, tmp_dy, tmp_ddy, tmap, *tddmap;
	int ix;
	double a,b,a1,b1,a2,b2;

	xm = CMAPDIM/2;
	p = CMAPDIM;

	d1y = 0.;
	d2y = 0.;
	d12y = 0.;

	memory->create(tmp_y,CMAPDIM*2,"cmap:tmp_y");
	memory->create(tmp_dy,CMAPDIM*2,"cmap:tmp_dy");
	memory->create(tmp_ddy,CMAPDIM*2,"cmap:tmp_ddy");
	memory->create(tmap,CMAPDIM2,CMAPDIM2,"cmap:tmap");
	memory->create(tddmap,CMAPDIM2,CMAPDIM2,"cmap:tddmap");

	// periodically expand the original map
	// use the expanded map for bicubic spline interpolation,
	// which is used to obtain the derivatives
	// actual interpolation is done with bicubic interpolation

	for (i = 0; i < CMAPDIM*2; i++) {
	ii = ((i+CMAPDIM-xm)%CMAPDIM);
	for (j = 0; j < CMAPDIM*2; j++) {
	jj = ((j+CMAPDIM-xm)%CMAPDIM);
	tmap[i][j] = map[ii][jj];
	}
	}

	for (i = 0; i < CMAPDIM*2; i++)
	spline(tmap[i], tddmap[i], CMAPDIM*2);

	for (i = xm; i < CMAPDIM+xm; i++) {
	phi = (i-xm)*CMAPDX-180.0;
	for (j = xm; j < CMAPDIM+xm; j++) {
	psi = (j-xm)*CMAPDX-180.0;
	ix = int((psi-CMAPXMIN)/CMAPDX);
	a = (CMAPXMIN+((ix+1)1.0)CMAPDX-psi)/CMAPDX;
	b = (psi-CMAPXMIN-((ix)1.0)CMAPDX)/CMAPDX;
	a1 = aaa-a;
	b1 = bbb-b;
	a2 = 3.0aa-1.0;
	b2 = 3.0bb-1.0;
	for (k = 0; k < CMAPDIM*2; k++) {
	tyyk = tmp_y[k];
	tdyk = tmp_dy[k];
	tyyk = atmap[k][ix]+btmap[k][ix+1]+
	(a1tddmap[k][ix]+b1tddmap[k][ix+1])(CMAPDXCMAPDX)/6.0;
	tdyk = (tmap[k][ix+1]-tmap[k][ix])/CMAPDX-
	(a2/6.0CMAPDXtddmap[k][ix])+(b2/6.0CMAPDXtddmap[k][ix+1]);
	tmp_y[k] = tyyk;
	tmp_dy[k] = tdyk;
	}

	spline(tmp_y,tmp_ddy,CMAPDIM+xm+xm);
	ix = int((phi-CMAPXMIN)/CMAPDX);
	a = (CMAPXMIN+((ix+1)1.0)CMAPDX-phi)/CMAPDX;
	b = (phi-CMAPXMIN-(ix1.0)CMAPDX)/CMAPDX;
	a1 = aaa-a;
	b1 = bbb-b;
	a2 = 3.0aa-1.0;
	b2 = 3.0bb-1.0;
	d1y = (tmp_y[ix+1]-tmp_y[ix])/CMAPDX-
	a2/6.0CMAPDXtmp_ddy[ix]+b2/6.0CMAPDXtmp_ddy[ix+1];
	spline(tmp_dy,tmp_ddy,CMAPDIM+xm+xm);
	ix = int((phi-CMAPXMIN)/CMAPDX);
	a = (CMAPXMIN+((ix+1)1.0)CMAPDX-phi)/CMAPDX;
	b = (phi-CMAPXMIN-(ix1.0)CMAPDX)/CMAPDX;
	a1 = aaa-a;
	b1 = bbb-b;
	a2 = 3.0aa-1.0;
	b2 = 3.0bb-1.0;
	d2y = atmp_dy[ix]+btmp_dy[ix+1]+
	(a1tmp_ddy[ix]+b1tmp_ddy[ix+1])(CMAPDXCMAPDX)/6.0;
	d12y = (tmp_dy[ix+1]-tmp_dy[ix])/CMAPDX-
	a2/6.0CMAPDXtmp_ddy[ix]+b2/6.0CMAPDXtmp_ddy[ix+1];
	d1yo[i%p][j%p] = d1y;
	d2yo[i%p][j%p] = d2y;
	d12yo[i%p][j%p] = d12y;
	}
	}

	memory->destroy(tmp_y);
	memory->destroy(tmp_dy);
	memory->destroy(tmp_ddy);
	memory->destroy(tmap);
	memory->destroy(tddmap);
	}

	/* ---------------------------------------------------------------------- */

	double FixCMAP::dihedral_angle_atan2(double fx, double fy, double fz,
	double ax, double ay, double az,
	double bx, double by, double bz,
	double absg)
	{
	// calculate the dihedral angle

	double angle, arg1, arg2;

	arg1 = absg(fxbx+fyby+fzbz);
	arg2 = axbx+ayby+az*bz;

	if (arg1 == 0 && arg2 == 0)
	error->all(FLERR,"CMAP: atan2 function cannot take 2 zero arguments");
	else {
	angle = atan2(arg1,arg2);
	angle = angle*180.0/MY_PI;
	}

	return angle;
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::bc_coeff(double gs, double d1gs, double d2gs, double d12gs)
	{
	// calculate the bicubic interpolation coefficients c_ij

	static int wt[16][16] =
	{ {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
	{-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0},
	{2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
	{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
	{0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1},
	{0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1},
	{-3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0},
	{9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2},
	{-6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2},
	{2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0},
	{-6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1},
	{4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1}
	};

	int i, j, k, in;
	double xx, x[16];

	for (i = 0; i < 4; i++) {
	x[i] = gs[i];
	x[i+4] = d1gs[i]*CMAPDX;
	x[i+8] = d2gs[i]*CMAPDX;
	x[i+12] = d12gs[i]CMAPDXCMAPDX;
	}

	in = 0;
	for (i = 0; i < 4; i++) {
	for (j = 0; j < 4; j++) {
	xx = 0.0;
	for (k = 0; k < 16; k++) xx += wt[in][k]*x[k];
	in++;
	cij[i][j] = xx;
	}
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixCMAP::bc_interpol(double x1, double x2, int low1, int low2, double *gs,
	double d1gs, double d2gs, double *d12gs)
	{
	// for a given point of interest and its corresponding grid square values,
	// gradients and cross-derivatives
	// calculate the interpolated value of the point of interest (POI)

	int i;
	double t, u, gs1l, gs2l;

	// set the interpolation coefficients

	bc_coeff(gs,d1gs,d2gs,d12gs);

	gs1l = g_axis[low1];
	gs2l = g_axis[low2];

	t = (x1-gs1l)/CMAPDX;
	u = (x2-gs2l)/CMAPDX;

	E = dEdPhi = dEdPsi = 0.0;

	for (i = 3; i >= 0; i--) {
	E = tE + ((cij[i][3]u+cij[i][2])u+cij[i][1])u+cij[i][0];
	dEdPhi = udEdPhi + (3.0cij[3][i]t+2.0cij[2][i])*t+cij[1][i];
	dEdPsi = tdEdPsi + (3.0cij[i][3]u+2.0cij[i][2])*u+cij[i][1];
	}

	dEdPhi *= (180.0/MY_PI/CMAPDX);
	dEdPsi *= (180.0/MY_PI/CMAPDX);
	}

	// ----------------------------------------------------------------------
	// ----------------------------------------------------------------------
	// methods to read and write data file
	// ----------------------------------------------------------------------
	// ----------------------------------------------------------------------

	void FixCMAP::read_data_header(char *line)
	{
	if (strstr(line,"crossterms")) {
	sscanf(line,BIGINT_FORMAT,&ncmap);
	} else error->all(FLERR,"Invalid read data header line for fix cmap");

	// didn't set in constructor b/c this fix could be defined
	// before newton command

	newton_bond = force->newton_bond;
	}

	/* ----------------------------------------------------------------------
	unpack N lines in buf from section of data file labeled by keyword
	id_offset is applied to atomID fields if multiple data files are read
	store CMAP interactions as if newton_bond = OFF, even if actually ON
	------------------------------------------------------------------------- */

	void FixCMAP::read_data_section(char keyword, int n, char buf,
	tagint id_offset)
	{
	int m,tmp,itype;
	tagint atom1,atom2,atom3,atom4,atom5;
	char *next;

	next = strchr(buf,'\n');
	*next = '\0';
	int nwords = atom->count_words(buf);
	*next = '\n';

	if (nwords != 7) {
	char str[128];
	sprintf(str,"Incorrect %s format in data file",keyword);
	error->all(FLERR,str);
	}

	// loop over lines of CMAP crossterms
	// tokenize the line into values
	// add crossterm to one of my atoms, depending on newton_bond

	for (int i = 0; i < n; i++) {
	next = strchr(buf,'\n');
	*next = '\0';
	sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
	" " TAGINT_FORMAT " " TAGINT_FORMAT,
	&tmp,&itype,&atom1,&atom2,&atom3,&atom4,&atom5);

	atom1 += id_offset;
	atom2 += id_offset;
	atom3 += id_offset;
	atom4 += id_offset;
	atom5 += id_offset;

	if ((m = atom->map(atom1)) >= 0) {
	if (num_crossterm[m] == CMAPMAX)
	error->one(FLERR,"Too many CMAP crossterms for one atom");
	crossterm_type[m][num_crossterm[m]] = itype;
	crossterm_atom1[m][num_crossterm[m]] = atom1;
	crossterm_atom2[m][num_crossterm[m]] = atom2;
	crossterm_atom3[m][num_crossterm[m]] = atom3;
	crossterm_atom4[m][num_crossterm[m]] = atom4;
	crossterm_atom5[m][num_crossterm[m]] = atom5;
	num_crossterm[m]++;
	}

	if ((m = atom->map(atom2)) >= 0) {
	if (num_crossterm[m] == CMAPMAX)
	error->one(FLERR,"Too many CMAP crossterms for one atom");
	crossterm_type[m][num_crossterm[m]] = itype;
	crossterm_atom1[m][num_crossterm[m]] = atom1;
	crossterm_atom2[m][num_crossterm[m]] = atom2;
	crossterm_atom3[m][num_crossterm[m]] = atom3;
	crossterm_atom4[m][num_crossterm[m]] = atom4;
	crossterm_atom5[m][num_crossterm[m]] = atom5;
	num_crossterm[m]++;
	}

	if ((m = atom->map(atom3)) >= 0) {
	if (num_crossterm[m] == CMAPMAX)
	error->one(FLERR,"Too many CMAP crossterms for one atom");
	crossterm_type[m][num_crossterm[m]] = itype;
	crossterm_atom1[m][num_crossterm[m]] = atom1;
	crossterm_atom2[m][num_crossterm[m]] = atom2;
	crossterm_atom3[m][num_crossterm[m]] = atom3;
	crossterm_atom4[m][num_crossterm[m]] = atom4;
	crossterm_atom5[m][num_crossterm[m]] = atom5;
	num_crossterm[m]++;
	}

	if ((m = atom->map(atom4)) >= 0) {
	if (num_crossterm[m] == CMAPMAX)
	error->one(FLERR,"Too many CMAP crossterms for one atom");
	crossterm_type[m][num_crossterm[m]] = itype;
	crossterm_atom1[m][num_crossterm[m]] = atom1;
	crossterm_atom2[m][num_crossterm[m]] = atom2;
	crossterm_atom3[m][num_crossterm[m]] = atom3;
	crossterm_atom4[m][num_crossterm[m]] = atom4;
	crossterm_atom5[m][num_crossterm[m]] = atom5;
	num_crossterm[m]++;
	}

	if ((m = atom->map(atom5)) >= 0) {
	if (num_crossterm[m] == CMAPMAX)
	error->one(FLERR,"Too many CMAP crossterms for one atom");
	crossterm_type[m][num_crossterm[m]] = itype;
	crossterm_atom1[m][num_crossterm[m]] = atom1;
	crossterm_atom2[m][num_crossterm[m]] = atom2;
	crossterm_atom3[m][num_crossterm[m]] = atom3;
	crossterm_atom4[m][num_crossterm[m]] = atom4;
	crossterm_atom5[m][num_crossterm[m]] = atom5;
	num_crossterm[m]++;
	}

	buf = next + 1;
	}
	}

	/* ---------------------------------------------------------------------- */

	bigint FixCMAP::read_data_skip_lines(char *keyword)
	{
	return ncmap;
	}

	/* ----------------------------------------------------------------------
	write Mth header line to file
	only called by proc 0
	------------------------------------------------------------------------- */

	void FixCMAP::write_data_header(FILE *fp, int mth)
	{
	fprintf(fp,BIGINT_FORMAT " cmap crossterms\n",ncmap);
	}

	/* ----------------------------------------------------------------------
	return size I own for Mth data section
	# of data sections = 1 for this fix
	nx = # of crossterms owned by my local atoms
	if newton_bond off, atom only owns crossterm if it is atom3
	ny = columns = type + 5 atom IDs
	------------------------------------------------------------------------- */

	void FixCMAP::write_data_section_size(int mth, int &nx, int &ny)
	{
	int i,m;

	tagint *tag = atom->tag;
	int nlocal = atom->nlocal;

	nx = 0;
	for (i = 0; i < nlocal; i++)
	for (m = 0; m < num_crossterm[i]; m++)
	if (crossterm_atom3[i][m] == tag[i]) nx++;

	ny = 6;
	}

	/* ----------------------------------------------------------------------
	pack values for Mth data section into 2d buf
	buf allocated by caller as owned crossterms by 6
	------------------------------------------------------------------------- */

	void FixCMAP::write_data_section_pack(int mth, double **buf)
	{
	int i,m;

	// 1st column = CMAP type
	// 2nd-6th columns = 5 atom IDs

	tagint *tag = atom->tag;
	int nlocal = atom->nlocal;

	int n = 0;
	for (i = 0; i < nlocal; i++) {
	for (m = 0; m < num_crossterm[i]; m++) {
	if (crossterm_atom3[i][m] != tag[i]) continue;
	buf[n][0] = ubuf(crossterm_type[i][m]).d;
	buf[n][1] = ubuf(crossterm_atom1[i][m]).d;
	buf[n][2] = ubuf(crossterm_atom2[i][m]).d;
	buf[n][3] = ubuf(crossterm_atom3[i][m]).d;
	buf[n][4] = ubuf(crossterm_atom4[i][m]).d;
	buf[n][5] = ubuf(crossterm_atom5[i][m]).d;
	n++;
	}
	}
	}

	/* ----------------------------------------------------------------------
	write section keyword for Mth data section to file
	use Molecules or Charges if that is only field, else use fix ID
	only called by proc 0
	------------------------------------------------------------------------- */

	void FixCMAP::write_data_section_keyword(int mth, FILE *fp)
	{
	fprintf(fp,"\nCMAP\n\n");
	}

	/* ----------------------------------------------------------------------
	write N lines from buf to file
	convert buf fields to int or double depending on styles
	index can be used to prepend global numbering
	only called by proc 0
	------------------------------------------------------------------------- */

	void FixCMAP::write_data_section(int mth, FILE *fp,
	int n, double **buf, int index)
	{
	for (int i = 0; i < n; i++)
	fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
	" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT "\n",
	index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
	(tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i,
	(tagint) ubuf(buf[i][4]).i,(tagint) ubuf(buf[i][5]).i);
	}

	// ----------------------------------------------------------------------
	// ----------------------------------------------------------------------
	// methods for restart and communication
	// ----------------------------------------------------------------------
	// ----------------------------------------------------------------------

	/* ----------------------------------------------------------------------
	pack entire state of Fix into one write
	------------------------------------------------------------------------- */

	void FixCMAP::write_restart(FILE *fp)
	{
	if (comm->me == 0) {
	int size = sizeof(bigint);
	fwrite(&size,sizeof(int),1,fp);
	fwrite(&ncmap,sizeof(bigint),1,fp);
	}
	}

	/* ----------------------------------------------------------------------
	use state info from restart file to restart the Fix
	------------------------------------------------------------------------- */

	void FixCMAP::restart(char *buf)
	{
	ncmap = ((bigint ) buf);
	}

	/* ----------------------------------------------------------------------
	pack values in local atom-based arrays for restart file
	------------------------------------------------------------------------- */

	int FixCMAP::pack_restart(int i, double *buf)
	{
	int n = 1;
	for (int m = 0; m < num_crossterm[i]; m++) {
	buf[n++] = ubuf(MAX(crossterm_type[i][m],-crossterm_type[i][m])).d;
	buf[n++] = ubuf(crossterm_atom1[i][m]).d;
	buf[n++] = ubuf(crossterm_atom2[i][m]).d;
	buf[n++] = ubuf(crossterm_atom3[i][m]).d;
	buf[n++] = ubuf(crossterm_atom4[i][m]).d;
	buf[n++] = ubuf(crossterm_atom5[i][m]).d;
	}
	buf[0] = n;

	return n;
	}

	/* ----------------------------------------------------------------------
	unpack values from atom->extra array to restart the fix
	------------------------------------------------------------------------- */

	void FixCMAP::unpack_restart(int nlocal, int nth)
	{
	double **extra = atom->extra;

	// skip to Nth set of extra values

	int n = 0;
	for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]);

	int count = static_cast<int> (extra[nlocal][n++]);
	num_crossterm[nlocal] = (count-1)/6;

	for (int m = 0; m < num_crossterm[nlocal]; m++) {
	crossterm_type[nlocal][m] = (int) ubuf(extra[nlocal][n++]).i;
	crossterm_atom1[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
	crossterm_atom2[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
	crossterm_atom3[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
	crossterm_atom4[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
	crossterm_atom5[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
	}
	}

	/* ----------------------------------------------------------------------
	maxsize of any atom's restart data
	------------------------------------------------------------------------- */

	int FixCMAP::maxsize_restart()
	{
	return 1 + CMAPMAX*6;
	}

	/* ----------------------------------------------------------------------
	size of atom nlocal's restart data
	------------------------------------------------------------------------- */

	int FixCMAP::size_restart(int nlocal)
	{
	return 1 + num_crossterm[nlocal]*6;
	}

	/* ----------------------------------------------------------------------
	allocate atom-based array
	------------------------------------------------------------------------- */

	void FixCMAP::grow_arrays(int nmax)
	{
	num_crossterm = memory->grow(num_crossterm,nmax,"cmap:num_crossterm");
	crossterm_type = memory->grow(crossterm_type,nmax,CMAPMAX,
	"cmap:crossterm_type");
	crossterm_atom1 = memory->grow(crossterm_atom1,nmax,CMAPMAX,
	"cmap:crossterm_atom1");
	crossterm_atom2 = memory->grow(crossterm_atom2,nmax,CMAPMAX,
	"cmap:crossterm_atom2");
	crossterm_atom3 = memory->grow(crossterm_atom3,nmax,CMAPMAX,
	"cmap:crossterm_atom3");
	crossterm_atom4 = memory->grow(crossterm_atom4,nmax,CMAPMAX,
	"cmap:crossterm_atom4");
	crossterm_atom5 = memory->grow(crossterm_atom5,nmax,CMAPMAX,
	"cmap:crossterm_atom5");

	// must initialize num_crossterm to 0 for added atoms
	// may never be set for some atoms when data file is read

	for (int i = nmax_previous; i < nmax; i++) num_crossterm[i] = 0;
	nmax_previous = nmax;
	}

	/* ----------------------------------------------------------------------
	copy values within local atom-based array
	------------------------------------------------------------------------- */

	void FixCMAP::copy_arrays(int i, int j, int delflag)
	{
	num_crossterm[j] = num_crossterm[i];

	for (int k = 0; k < num_crossterm[j]; k++){
	crossterm_type[j][k] = crossterm_type[i][k];
	crossterm_atom1[j][k] = crossterm_atom1[i][k];
	crossterm_atom2[j][k] = crossterm_atom2[i][k];
	crossterm_atom3[j][k] = crossterm_atom3[i][k];
	crossterm_atom4[j][k] = crossterm_atom4[i][k];
	crossterm_atom5[j][k] = crossterm_atom5[i][k];
	}
	}

	/* ----------------------------------------------------------------------
	initialize one atom's array values, called when atom is created
	------------------------------------------------------------------------- */

	void FixCMAP::set_arrays(int i)
	{
	num_crossterm[i] = 0;
	}

	/* ----------------------------------------------------------------------
	pack values in local atom-based array for exchange with another proc
	------------------------------------------------------------------------- */

	int FixCMAP::pack_exchange(int i, double *buf)
	{
	int n = 0;
	buf[n++] = ubuf(num_crossterm[i]).d;
	for (int m = 0; m < num_crossterm[i]; m++) {
	buf[n++] = ubuf(crossterm_type[i][m]).d;
	buf[n++] = ubuf(crossterm_atom1[i][m]).d;
	buf[n++] = ubuf(crossterm_atom2[i][m]).d;
	buf[n++] = ubuf(crossterm_atom3[i][m]).d;
	buf[n++] = ubuf(crossterm_atom4[i][m]).d;
	buf[n++] = ubuf(crossterm_atom5[i][m]).d;
	}
	return n;
	}

	/* ----------------------------------------------------------------------
	unpack values in local atom-based array from exchange with another proc
	------------------------------------------------------------------------- */

	int FixCMAP::unpack_exchange(int nlocal, double *buf)
	{
	int n = 0;
	num_crossterm[nlocal] = (int) ubuf(buf[n++]).i;
	for (int m = 0; m < num_crossterm[nlocal]; m++) {
	crossterm_type[nlocal][m] = (int) ubuf(buf[n++]).i;
	crossterm_atom1[nlocal][m] = (tagint) ubuf(buf[n++]).i;
	crossterm_atom2[nlocal][m] = (tagint) ubuf(buf[n++]).i;
	crossterm_atom3[nlocal][m] = (tagint) ubuf(buf[n++]).i;
	crossterm_atom4[nlocal][m] = (tagint) ubuf(buf[n++]).i;
	crossterm_atom5[nlocal][m] = (tagint) ubuf(buf[n++]).i;
	}
	return n;
	}

	/* ----------------------------------------------------------------------
	memory usage of local atom-based arrays
	------------------------------------------------------------------------- */

	double FixCMAP::memory_usage()
	{
	int nmax = atom->nmax;
	double bytes = nmax * sizeof(int); // num_crossterm
	bytes += nmaxCMAPMAX sizeof(int); // crossterm_type
	bytes += 5nmaxCMAPMAX * sizeof(int); // crossterm_atom 12345
	bytes += maxcrossterm6 sizeof(int); // crosstermlist
	return bytes;
	}

fix_cmap.cppNo OneTemporaryActions

File Metadata

fix_cmap.cppView Options

Event Timeline

fix_cmap.cpp
No OneTemporary
Actions

fix_cmap.cpp
View Options