Diffusion lammps (lammm-master-deprecated) Browse

dihedral_table.cpp
56933a0bd14alammm-master-deprecated

dihedral_table.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 author: Andrew Jewett (jewett.aij g m ail)
			The cyclic tridiagonal matrix solver was borrowed from
			the "tridiag.c" written by Gerard Jungman for GSL
			------------------------------------------------------------------------- */

			#include <cmath>
			#include <cstdlib>
			#include <cstring>
			#include <string>
			#include <fstream>
			#include <iostream>
			#define NDEBUG //(disable assert())
			#include <cassert>
			#include <sstream>
			using namespace std;

			#include "atom.h"
			#include "comm.h"
			#include "neighbor.h"
			#include "domain.h"
			#include "force.h"
			#include "update.h"
			#include "memory.h"
			#include "error.h"
			#include "dihedral_table.h"

			// Additional header files available for numerical debugging:
			//#undef NDEBUG
			//#define DIH_DEBUG_NUM
			//#ifdef DIH_DEBUG_NUM
			//#include "dihedral_table_dbg.h" //in "supporting_files/debug_numerical/"
			//#endif
			// and pointless posturing:
			//#include "dihedral_table_nd_mod.h" //in "supporting_files/nd/"

			using namespace LAMMPS_NS;
			using namespace DIHEDRAL_TABLE_NS;

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

			DihedralTable::DihedralTable(LAMMPS *lmp) : Dihedral(lmp)
			{
			ntables = 0;
			tables = NULL;
			}

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

			DihedralTable::~DihedralTable()
			{
			for (int m = 0; m < ntables; m++) free_table(&tables[m]);
			memory->sfree(tables);

			if (allocated) {
			memory->destroy(setflag);
			//memory->destroy(phi0); <- equilibrium angles not supported
			memory->destroy(tabindex);
			}
			}

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

			void DihedralTable::compute(int eflag, int vflag)
			{
			int i1,i2,i3,i4,n,type;
			double edihedral,f1[3],f2[3],f3[3],f4[3];

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

			int **dihedrallist = neighbor->dihedrallist;
			int ndihedrallist = neighbor->ndihedrallist;
			int nlocal = atom->nlocal;
			int newton_bond = force->newton_bond;

			// The dihedral angle "phi" is the angle between n123 and n234
			// the planes defined by atoms i1,i2,i3, and i2,i3,i4.
			//
			// Definitions of vectors: vb12, vb23, vb34, perp12on23
			// proj12on23, perp43on32, proj43on32
			//
			// Note: The positions of the 4 atoms are labeled x[i1], x[i2], x[i3], x[i4]
			// (which are also vectors)
			//
			// proj12on23 proj34on23
			// ---------> ----------->
			// .
			// .
			// .
			// x[i2] . x[i3]
			// . __@----------vb23-------->@ . . . . .
			// /\|\ /\| \ \|
			// \| / \ \|
			// \| / \ \|
			// perp12vs23 / \ \|
			// \| / \ perp34vs23
			// \| vb12 \ \|
			// \| / vb34 \|
			// \| / \ \|
			// \| / \ \|
			// \| / \ \|
			// @ \ \|
			// _\\| \\|/
			// x[i1] @
			//
			// x[i4]
			//

			double vb12[g_dim]; // displacement vector from atom i1 towards atom i2
			// vb12[d] = x[i2][d] - x[i1][d] (for d=0,1,2)
			double vb23[g_dim]; // displacement vector from atom i2 towards atom i3
			// vb23[d] = x[i3][d] - x[i2][d] (for d=0,1,2)
			double vb34[g_dim]; // displacement vector from atom i3 towards atom i4
			// vb34[d] = x[i4][d] - x[i3][d] (for d=0,1,2)

			// n123 & n234: These two unit vectors are normal to the planes
			// defined by atoms 1,2,3 and 2,3,4.
			double n123[g_dim]; //n123=vb23 x vb12 / \|vb23 x vb12\| ("x" is cross product)
			double n234[g_dim]; //n234=vb23 x vb34 / \|vb23 x vb34\| ("x" is cross product)

			double proj12on23[g_dim];
			// proj12on23[d] = (vb23[d]/\|vb23\|) * DotProduct(vb12,vb23)/\|vb12\|*\|vb23\|
			double proj34on23[g_dim];
			// proj34on23[d] = (vb34[d]/\|vb23\|) * DotProduct(vb34,vb23)/\|vb34\|*\|vb23\|
			double perp12on23[g_dim];
			// perp12on23[d] = v12[d] - proj12on23[d]
			double perp34on23[g_dim];
			// perp34on23[d] = v34[d] - proj34on23[d]


			edihedral = 0.0;
			if (eflag \|\| vflag) ev_setup(eflag,vflag);
			else evflag = 0;


			for (n = 0; n < ndihedrallist; n++) {

			i1 = dihedrallist[n][0];
			i2 = dihedrallist[n][1];
			i3 = dihedrallist[n][2];
			i4 = dihedrallist[n][3];
			type = dihedrallist[n][4];

			// ------ Step 1: Compute the dihedral angle "phi" ------
			//

			// Phi() calculates the dihedral angle.
			// This function also calculates the vectors:
			// vb12, vb23, vb34, n123, and n234, which we will need later.


			double phi = Phi(x[i1], x[i2], x[i3], x[i4], domain,
			vb12, vb23, vb34, n123, n234);


			// ------ Step 2: Compute the gradient of phi with atomic position: ------
			//
			// Gradient variables:
			//
			// dphi_dx1, dphi_dx2, dphi_dx3, dphi_dx4 are the gradients of phi with
			// respect to the atomic positions of atoms i1, i2, i3, i4, respectively.
			// As an example, consider dphi_dx1. The d'th element is:
			double dphi_dx1[g_dim]; // d phi
			double dphi_dx2[g_dim]; // dphi_dx1[d] = ---------- (partial derivatives)
			double dphi_dx3[g_dim]; // d x[i1][d]
			double dphi_dx4[g_dim]; //where d=0,1,2 corresponds to x,y,z (if g_dim==3)

			double dot123 = DotProduct(vb12, vb23);
			double dot234 = DotProduct(vb23, vb34);
			double L23sqr = DotProduct(vb23, vb23);
			double L23 = sqrt(L23sqr); // (central bond length)
			double inv_L23sqr = 0.0;
			double inv_L23 = 0.0;
			if (L23sqr != 0.0) {
			inv_L23sqr = 1.0 / L23sqr;
			inv_L23 = 1.0 / L23;
			}
			double neg_inv_L23 = -inv_L23;
			double dot123_over_L23sqr = dot123 * inv_L23sqr;
			double dot234_over_L23sqr = dot234 * inv_L23sqr;

			for (int d=0; d < g_dim; ++d) {
			// See figure above for a visual definitions of these vectors:
			proj12on23[d] = vb23[d] * dot123_over_L23sqr;
			proj34on23[d] = vb23[d] * dot234_over_L23sqr;
			perp12on23[d] = vb12[d] - proj12on23[d];
			perp34on23[d] = vb34[d] - proj34on23[d];
			}


			// --- Compute the gradient vectors dphi/dx1 and dphi/dx4: ---

			// These two gradients point in the direction of n123 and n234,
			// and are scaled by the distances of atoms 1 and 4 from the central axis.
			// Distance of atom 1 to central axis:
			double perp12on23_len = sqrt(DotProduct(perp12on23, perp12on23));
			// Distance of atom 4 to central axis:
			double perp34on23_len = sqrt(DotProduct(perp34on23, perp34on23));

			double inv_perp12on23 = 0.0;
			if (perp12on23_len != 0.0) inv_perp12on23 = 1.0 / perp12on23_len;
			double inv_perp34on23 = 0.0;
			if (perp34on23_len != 0.0) inv_perp34on23 = 1.0 / perp34on23_len;

			for (int d=0; d < g_dim; ++d) {
			dphi_dx1[d] = n123[d] * inv_perp12on23;
			dphi_dx4[d] = n234[d] * inv_perp34on23;
			}

			// --- Compute the gradient vectors dphi/dx2 and dphi/dx3: ---
			//
			// This is more tricky because atoms 2 and 3 are shared by both planes
			// 123 and 234 (the angle between which defines "phi"). Moving either
			// one of these atoms effects both the 123 and 234 planes
			// Both the 123 and 234 planes intersect with the plane perpendicular to the
			// central bond axis (vb23). The two lines where these intersections occur
			// will shift when you move either atom 2 or atom 3. The angle between
			// these lines is the dihedral angle, phi. We can define four quantities:
			// dphi123_dx2 is the change in "phi" due to the movement of the 123 plane
			// ...as a result of moving atom 2.
			// dphi234_dx2 is the change in "phi" due to the movement of the 234 plane
			// ...as a result of moving atom 2.
			// dphi123_dx3 is the change in "phi" due to the movement of the 123 plane
			// ...as a result of moving atom 3.
			// dphi234_dx3 is the change in "phi" due to the movement of the 234 plane
			// ...as a result of moving atom 3.

			double proj12on23_len = dot123 * inv_L23;
			double proj34on23_len = dot234 * inv_L23;
			// Interpretation:
			//The magnitude of "proj12on23_len" is the length of the proj12on23 vector.
			//The sign is positive if it points in the same direction as the central
			//bond (vb23). Otherwise it is negative. The same goes for "proj34on23".
			//(In the example figure in the comment above, both variables are positive.)

			// The forumula used in the 8 lines below explained here:
			// "supporting_information/doc/gradient_formula_explanation/"
			double dphi123_dx2_coef = neg_inv_L23 * (L23 + proj12on23_len);
			double dphi234_dx2_coef = inv_L23 * proj34on23_len;

			double dphi234_dx3_coef = neg_inv_L23 * (L23 + proj34on23_len);
			double dphi123_dx3_coef = inv_L23 * proj12on23_len;

			for (int d=0; d < g_dim; ++d) {
			// Recall that the n123 and n234 plane normal vectors are proportional to
			// the dphi/dx1 and dphi/dx2 gradients vectors
			// It turns out we can save slightly more CPU cycles by expressing
			// dphi/dx2 and dphi/dx3 as linear combinations of dphi/dx1 and dphi/dx2
			// which we computed already (instead of n123 & n234).
			dphi_dx2[d] = dphi123_dx2_coefdphi_dx1[d] + dphi234_dx2_coefdphi_dx4[d];
			dphi_dx3[d] = dphi123_dx3_coefdphi_dx1[d] + dphi234_dx3_coefdphi_dx4[d];
			}




			#ifdef DIH_DEBUG_NUM
			// ----- Numerical test? -----

			cerr << " -- testing gradient for dihedral (n="<<n<<") for atoms ("
			<< i1 << "," << i2 << "," << i3 << "," << i4 << ") --" << endl;

			PrintGradientComparison(*this, dphi_dx1, dphi_dx2, dphi_dx3, dphi_dx4,
			domain, x[i1], x[i2], x[i3], x[i4]);

			for (int d=0; d < g_dim; ++d) {
			// The sum of all the gradients should be near 0. (translational symmetry)
			cerr <<"sum_gradients["<<d<<"]="<<dphi_dx1[d]<<"+"<<dphi_dx2[d]<<"+"<<dphi_dx3[d]<<"+"<<dphi_dx4[d]<<"="<<dphi_dx1[d]+dphi_dx2[d]+dphi_dx3[d]+dphi_dx4[d]<<endl;
			// These should sum to zero
			assert(abs(dphi_dx1[d]+dphi_dx2[d]+dphi_dx3[d]+dphi_dx4[d]) < 0.0002/L23);
			}
			#endif // #ifdef DIH_DEBUG_NUM




			// ----- Step 3: Calculate the energy and force in the phi direction -----

			// tabulated force & energy
			double u, m_du_dphi; //u = energy. m_du_dphi = "minus" du/dphi
			assert((0.0 <= phi) && (phi <= TWOPI));

			uf_lookup(type, phi, u, m_du_dphi);


			if (eflag) edihedral = u;

			// ----- Step 4: Calculate the force direction in real space -----

			// chain rule:
			// d U d U d phi
			// -f = ----- = ----- * -----
			// d x d phi d x
			for(int d=0; d < g_dim; ++d) {
			f1[d] = m_du_dphi * dphi_dx1[d];
			f2[d] = m_du_dphi * dphi_dx2[d];
			f3[d] = m_du_dphi * dphi_dx3[d];
			f4[d] = m_du_dphi * dphi_dx4[d];
			}

			// apply force to each of 4 atoms

			if (newton_bond \|\| i1 < nlocal) {
			f[i1][0] += f1[0];
			f[i1][1] += f1[1];
			f[i1][2] += f1[2];
			}

			if (newton_bond \|\| i2 < nlocal) {
			f[i2][0] += f2[0];
			f[i2][1] += f2[1];
			f[i2][2] += f2[2];
			}

			if (newton_bond \|\| i3 < nlocal) {
			f[i3][0] += f3[0];
			f[i3][1] += f3[1];
			f[i3][2] += f3[2];
			}

			if (newton_bond \|\| i4 < nlocal) {
			f[i4][0] += f4[0];
			f[i4][1] += f4[1];
			f[i4][2] += f4[2];
			}

			if (evflag)
			ev_tally(i1,i2,i3,i4,
			nlocal,newton_bond,edihedral,
			f1,f3,f4,
			vb12[0],vb12[1],vb12[2],
			vb23[0],vb23[1],vb23[2],
			vb34[0],vb34[1],vb34[2]);
			}
			} // void DihedralTable::compute()







			// single() calculates the dihedral-angle energy of atoms i1, i2, i3, i4.
			double DihedralTable::single(int type, int i1, int i2, int i3, int i4)
			{
			double vb12[g_dim];
			double vb23[g_dim];
			double vb34[g_dim];
			double n123[g_dim];
			double n234[g_dim];

			double **x = atom->x;

			double phi = Phi(x[i1], x[i2], x[i3], x[i4], domain,
			vb12, vb23, vb34, n123, n234);

			double u;
			u_lookup(type, phi, u); //Calculate the energy, and store it in "u"

			return u;
			}


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



			void DihedralTable::allocate()
			{
			allocated = 1;
			int n = atom->ndihedraltypes;

			memory->create(tabindex,n+1,"dihedral:tabindex");
			//memory->create(phi0,n+1,"dihedral:phi0"); <-equilibrium angles not supported
			memory->create(setflag,n+1,"dihedral:setflag");
			for (int i = 1; i <= n; i++) setflag[i] = 0;
			}


			/* ----------------------------------------------------------------------
			global settings
			------------------------------------------------------------------------- */

			void DihedralTable::settings(int narg, char **arg)
			{
			if (narg != 2) error->all(FLERR,"Illegal dihedral_style command");

			if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
			else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
			else error->all(FLERR,"Unknown table style in dihedral style table");

			tablength = force->inumeric(arg[1]);
			if (tablength < 3)
			error->all(FLERR,"Illegal number of dihedral table entries");
			// delete old tables, since cannot just change settings

			for (int m = 0; m < ntables; m++) free_table(&tables[m]);
			memory->sfree(tables);

			if (allocated) {
			memory->destroy(setflag);
			memory->destroy(tabindex);
			}
			allocated = 0;

			ntables = 0;
			tables = NULL;
			}



			/* ----------------------------------------------------------------------
			set coeffs for one type
			------------------------------------------------------------------------- */


			void DihedralTable::coeff(int narg, char **arg)
			{
			if (narg != 3) error->all(FLERR,"Illegal dihedral_coeff command");
			if (!allocated) allocate();

			int ilo,ihi;
			force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi);

			int me;
			MPI_Comm_rank(world,&me);
			tables = (Table *)
			memory->srealloc(tables,(ntables+1)*sizeof(Table), "dihedral:tables");
			Table *tb = &tables[ntables];
			null_table(tb);
			if (me == 0) read_table(tb,arg[1],arg[2]);
			bcast_table(tb);


			// --- check the angle data for range errors ---
			// --- and resolve issues with periodicity ---

			if (tb->ninput < 3) {
			string err_msg;
			err_msg = string("Invalid dihedral table length (")
			+ string(arg[2]) + string(").");
			error->one(FLERR,err_msg.c_str());
			}

			// check for monotonicity
			for (int i=0; i < tb->ninput-1; i++) {
			if (tb->phifile[i] >= tb->phifile[i+1]) {
			stringstream i_str;
			i_str << i+1;
			string err_msg =
			string("Dihedral table values are not increasing (") +
			string(arg[2]) + string(", ")+i_str.str()+string("th entry)");
			if (i==0)
			err_msg += string("\n(This is probably a mistake with your table format.)\n");
			error->all(FLERR,err_msg.c_str());
			}
			}

			// check the range of angles
			double philo = tb->phifile[0];
			double phihi = tb->phifile[tb->ninput-1];
			if (tb->use_degrees) {
			if ((phihi - philo) >= 360) {
			string err_msg;
			err_msg = string("Dihedral table angle range must be < 360 degrees (")
			+string(arg[2]) + string(").");
			error->all(FLERR,err_msg.c_str());
			}
			}
			else {
			if ((phihi - philo) >= TWOPI) {
			string err_msg;
			err_msg = string("Dihedral table angle range must be < 2*PI radians (")
			+ string(arg[2]) + string(").");
			error->all(FLERR,err_msg.c_str());
			}
			}

			// convert phi from degrees to radians
			if (tb->use_degrees) {
			for (int i=0; i < tb->ninput; i++) {
			tb->phifile[i] *= PI/180.0;
			// I assume that if angles are in degrees, then the forces (f=dU/dphi)
			// are specified with "phi" in radians as well.
			tb->ffile[i] *= 180.0/PI;
			}
			}

			// We want all the phi dihedral angles to lie in the range from 0 to 2*PI.
			// But I don't want to restrict users to input their data in this range.
			// We also want the angles to be sorted in increasing order.
			// This messy code fixes these problems with the user's data:
			{
			double *phifile_tmp = new double [tb->ninput]; //temporary arrays
			double *ffile_tmp = new double [tb->ninput]; //used for sorting
			double *efile_tmp = new double [tb->ninput];

			// After re-imaging, does the range of angles cross the 0 or 2*PI boundary?
			// If so, find the discontinuity:
			int i_discontinuity = tb->ninput;
			for (int i=0; i < tb->ninput; i++) {
			double phi = tb->phifile[i];
			// Add a multiple of 2PI to phi until it lies in the range [0, 2PI).
			phi -= TWOPI * floor(phi/TWOPI);
			phifile_tmp[i] = phi;
			efile_tmp[i] = tb->efile[i];
			ffile_tmp[i] = tb->ffile[i];
			if ((i>0) && (phifile_tmp[i] < phifile_tmp[i-1])) {
			//There should only be at most one discontinuity, because we have
			//insured that the data was sorted before imaging, and because the
			//range of angle values does not exceed 2*PI.
			assert(i_discontinuity == tb->ninput); //<-should only happen once
			i_discontinuity = i;
			}
			}

			int I = 0;
			for (int i = i_discontinuity; i < tb->ninput; i++) {
			tb->phifile[I] = phifile_tmp[i];
			tb->efile[I] = efile_tmp[i];
			tb->ffile[I] = ffile_tmp[i];
			assert((I==0) \|\| (tb->phifile[I] > tb->phifile[I-1]));
			I++;
			}
			for (int i = 0; i < i_discontinuity; i++) {
			tb->phifile[I] = phifile_tmp[i];
			tb->efile[I] = efile_tmp[i];
			tb->ffile[I] = ffile_tmp[i];
			assert((I==0) \|\| (tb->phifile[I] > tb->phifile[I-1]));
			I++;
			}
			assert(I == tb->ninput);
			}

			// spline read-in and compute r,e,f vectors within table

			spline_table(tb);
			compute_table(tb);

			// Optional: allow the user to print out the interpolated spline tables

			if (me == 0) {
			if (checkU_fname && (strlen(checkU_fname) != 0))
			{
			ofstream checkU_file;
			checkU_file.open(checkU_fname, ios::out);
			for (int i=0; i < tablength; i++) {
			double phi = i*TWOPI/tablength;
			double u = tb->e[i];
			if (tb->use_degrees)
			phi *= 180.0/PI;
			checkU_file << phi << " " << u << "\n";
			}
			checkU_file.close();
			}
			if (checkF_fname && (strlen(checkF_fname) != 0))
			{
			ofstream checkF_file;
			checkF_file.open(checkF_fname, ios::out);
			for (int i=0; i < tablength; i++)
			{
			double phi = i*TWOPI/tablength;
			double f;
			if ((tabstyle == SPLINE) && (tb->f_unspecified)) {
			double dU_dphi =
			// (If the user did not specify the forces now, AND the user
			// selected the "spline" option, (as opposed to "linear")
			// THEN the tb->f array is uninitialized, so there's
			// no point to print out the contents of the tb->f[] array.
			// Instead, later on, we will calculate the force using the
			// -cyc_splintD() routine to calculate the derivative of the
			// energy spline, using the energy data (tb->e[]).
			// To be nice and report something, I do the same thing here.)
			cyc_splintD(tb->phi, tb->e, tb->e2, tablength, TWOPI,phi);
			f = -dU_dphi;
			}
			else
			// Otherwise we calculated the tb->f[] array. Report its contents.
			f = tb->f[i];
			if (tb->use_degrees) {
			phi *= 180.0/PI;
			// If the user wants degree angle units, we should convert our
			// internal force tables (in energy/radians) to (energy/degrees)
			f *= PI/180.0;
			}
			checkF_file << phi << " " << f << "\n";
			}
			checkF_file.close();
			} // if (checkF_fname && (strlen(checkF_fname) != 0))
			} // if (me == 0)

			// store ptr to table in tabindex
			int count = 0;
			for (int i = ilo; i <= ihi; i++)
			{
			tabindex[i] = ntables;
			//phi0[i] = tb->phi0; <- equilibrium dihedral angles not supported
			setflag[i] = 1;
			count++;
			}
			ntables++;

			if (count == 0)
			error->all(FLERR,"Illegal dihedral_coeff command");

			} //DihedralTable::coeff()


			/* ----------------------------------------------------------------------
			proc 0 writes to restart file
			------------------------------------------------------------------------- */

			void DihedralTable::write_restart(FILE *fp)
			{
			fwrite(&tabstyle,sizeof(int),1,fp);
			fwrite(&tablength,sizeof(int),1,fp);
			}

			/* ----------------------------------------------------------------------
			proc 0 reads from restart file, bcasts
			------------------------------------------------------------------------- */

			void DihedralTable::read_restart(FILE *fp)
			{
			if (comm->me == 0) {
			fread(&tabstyle,sizeof(int),1,fp);
			fread(&tablength,sizeof(int),1,fp);
			}

			//MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world); <-looks like a bug. Andrew 2012
			MPI_Bcast(&tabstyle,1,MPI_INT,0,world);
			MPI_Bcast(&tablength,1,MPI_INT,0,world);

			allocate();
			}


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

			void DihedralTable::null_table(Table *tb)
			{
			tb->phifile = tb->efile = tb->ffile = NULL;
			tb->e2file = tb->f2file = NULL;
			tb->phi = tb->e = tb->de = NULL;
			tb->f = tb->df = tb->e2 = tb->f2 = NULL;
			}

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

			void DihedralTable::free_table(Table *tb)
			{
			memory->destroy(tb->phifile);
			memory->destroy(tb->efile);
			memory->destroy(tb->ffile);
			memory->destroy(tb->e2file);
			memory->destroy(tb->f2file);

			memory->destroy(tb->phi);
			memory->destroy(tb->e);
			memory->destroy(tb->de);
			memory->destroy(tb->f);
			memory->destroy(tb->df);
			memory->destroy(tb->e2);
			memory->destroy(tb->f2);
			}

			/* ----------------------------------------------------------------------
			read table file, only called by proc 0
			------------------------------------------------------------------------- */

			void DihedralTable::read_table(Table tb, char file, char *keyword)
			{
			char line[MAXLINE];

			// open file

			FILE *fp = fopen(file,"r");
			if (fp == NULL) {
			string err_msg = string("Cannot open file ") + string(file);
			error->one(FLERR,err_msg.c_str());
			}

			// loop until section found with matching keyword

			while (1) {
			if (fgets(line,MAXLINE,fp) == NULL) {
			string err_msg=string("Did not find keyword \"")
			+string(keyword)+string("\" in dihedral table file.");
			error->one(FLERR, err_msg.c_str());
			}
			if (strspn(line," \t\n\r") == strlen(line)) continue; // blank line
			if (line[0] == '#') continue; // comment
			if (strstr(line,keyword) == line) break; // matching keyword
			fgets(line,MAXLINE,fp); // no match, skip section
			param_extract(tb,line);
			fgets(line,MAXLINE,fp);
			for (int i = 0; i < tb->ninput; i++)
			fgets(line,MAXLINE,fp);
			}

			// read args on 2nd line of section
			// allocate table arrays for file values

			fgets(line,MAXLINE,fp);
			param_extract(tb,line);
			memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
			memory->create(tb->efile,tb->ninput,"dihedral:efile");
			memory->create(tb->ffile,tb->ninput,"dihedral:ffile");

			// read a,e,f table values from file

			int itmp;
			for (int i = 0; i < tb->ninput; i++) {
			fgets(line,MAXLINE,fp);
			// Skip blank lines and delete text following a '#' character
			char *pe = strchr(line, '#');
			if (pe != NULL) *pe = '\0'; //terminate string at '#' character
			char *pc = line;
			while ((pc != '\0') && isspace(pc))
			pc++;
			if (*pc != '\0') { //If line is not a blank line
			stringstream line_ss(line);
			if (tb->f_unspecified) {
			//sscanf(line,"%d %lg %lg",
			// &itmp,&tb->phifile[i],&tb->efile[i]);
			line_ss >> itmp;
			line_ss >> tb->phifile[i];
			line_ss >> tb->efile[i];
			}
			else {
			//sscanf(line,"%d %lg %lg %lg",
			// &itmp,&tb->phifile[i],&tb->efile[i],&tb->ffile[i]);
			line_ss >> itmp;
			line_ss >> tb->phifile[i];
			line_ss >> tb->efile[i];
			line_ss >> tb->ffile[i];
			}
			if (! line_ss) {
			stringstream err_msg;
			err_msg << "Read error in table "<< keyword<<", near line "<<i+1<<"\n"
			<< " (Check to make sure the number of colums is correct.)";
			if ((! tb->f_unspecified) && (i==0))
			err_msg << "\n (This sometimes occurs if users forget to specify the \"NOF\" option.)\n";
			error->one(FLERR, err_msg.str().c_str());
			}
			}
			else //if it is a blank line, then skip it.
			i--;
			}

			fclose(fp);
			}

			/* ----------------------------------------------------------------------
			build spline representation of e,f over entire range of read-in table
			this function sets these values in e2file,f2file.
			I also perform a crude check for force & energy consistency.
			------------------------------------------------------------------------- */

			void DihedralTable::spline_table(Table *tb)
			{
			memory->create(tb->e2file,tb->ninput,"dihedral:e2file");
			memory->create(tb->f2file,tb->ninput,"dihedral:f2file");

			cyc_spline(tb->phifile, tb->efile, tb->ninput, TWOPI, tb->e2file);

			if (! tb->f_unspecified) {
			cyc_spline(tb->phifile, tb->ffile, tb->ninput, TWOPI, tb->f2file);
			}

			// CHECK to help make sure the user calculated forces in a way
			// which is grossly numerically consistent with the energy table.
			// --------------------------------------------------
			// This is an ugly piece of code
			// --------------------------------------------------
			if (! tb->f_unspecified) {
			int num_disagreements = 0;
			for (int i=0; i<tb->ninput; i++) {

			// Calculate what the force should be at the control points
			// by using linear interpolation of the derivatives of the energy:

			double phi_i = tb->phifile[i];

			// First deal with periodicity (I hate this part)
			double phi_im1, phi_ip1;
			int im1 = i-1;
			if (im1 < 0) {
			im1 += tb->ninput;
			phi_im1 = tb->phifile[im1] - TWOPI;
			}
			else
			phi_im1 = tb->phifile[im1];
			int ip1 = i+1;
			if (ip1 >= tb->ninput) {
			ip1 -= tb->ninput;
			phi_ip1 = tb->phifile[ip1] + TWOPI;
			}
			else
			phi_ip1 = tb->phifile[ip1];

			// Now calculate the midpoints above and below phi_i = tb->phifile[i]
			double phi_lo= 0.5*(phi_im1 + phi_i); //midpoint between phi_im1 and phi_i
			double phi_hi= 0.5*(phi_i + phi_ip1); //midpoint between phi_i and phi_ip1

			// Use a linear approximation to the derivative at these two midpoints
			double dU_dphi_lo =
			(tb->efile[i] - tb->efile[im1])
			/
			(phi_i - phi_im1);
			double dU_dphi_hi =
			(tb->efile[ip1] - tb->efile[i])
			/
			(phi_ip1 - phi_i);

			// Now calculate the derivative at position
			// phi_i (=tb->phifile[i]) using linear interpolation

			double a = (phi_i - phi_lo) / (phi_hi - phi_lo);
			double b = (phi_hi - phi_i) / (phi_hi - phi_lo);
			double dU_dphi = adU_dphi_lo + bdU_dphi_hi;
			double f = -dU_dphi;
			// Alternately, we could use spline interpolation instead:
			// double f = - splintD(tb->phifile, tb->efile, tb->e2file,
			// tb->ninput, TWOPI, tb->phifile[i]);
			// This is the way I originally did it, but I trust
			// the ugly simple linear way above better.
			// Recall this entire block of code doess not calculate
			// anything important. It does not have to be perfect.
			// We are only checking for stupid user errors here.

			if ((f != 0.0) &&
			(tb->ffile[i] != 0.0) &&
			((f/tb->ffile[i] < 0.5) \|\| (f/tb->ffile[i] > 2.0))) {
			num_disagreements++;
			}
			} // for (int i=0; i<tb->ninput; i++)

			if (num_disagreements > tb->ninput/2) {
			string msg("Dihedral table has inconsistent forces and energies. (Try \"NOF\".)\n");
			error->all(FLERR,msg.c_str());
			}

			} // check for consistency if (! tb->f_unspecified)

			} // DihedralTable::spline_table()


			/* ----------------------------------------------------------------------
			compute a,e,f vectors from splined values
			------------------------------------------------------------------------- */

			void DihedralTable::compute_table(Table *tb)
			{
			//delta = table spacing in dihedral angle for tablength (cyclic) bins
			tb->delta = TWOPI / tablength;
			tb->invdelta = 1.0/tb->delta;
			tb->deltasq6 = tb->delta*tb->delta / 6.0;

			// N evenly spaced bins in dihedral angle from 0 to 2*PI
			// phi,e,f = value at lower edge of bin
			// de,df values = delta values of e,f (cyclic, in this case)
			// phi,e,f,de,df are arrays containing "tablength" number of entries

			memory->create(tb->phi,tablength,"dihedral:phi");
			memory->create(tb->e,tablength,"dihedral:e");
			memory->create(tb->de,tablength,"dihedral:de");
			memory->create(tb->f,tablength,"dihedral:f");
			memory->create(tb->df,tablength,"dihedral:df");
			memory->create(tb->e2,tablength,"dihedral:e2");
			memory->create(tb->f2,tablength,"dihedral:f2");

			// Use cubic spline interpolation to calculate the entries in the
			// internal table. (This is true regardless...even if tabstyle!=SPLINE.)
			for (int i = 0; i < tablength; i++) {
			double phi = i*tb->delta;
			tb->phi[i] = phi;
			tb->e[i]= cyc_splint(tb->phifile,tb->efile,tb->e2file,tb->ninput,TWOPI,phi);
			if (! tb->f_unspecified)
			tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,TWOPI,phi);
			}

			if (tabstyle == LINEAR) {
			if (tb->f_unspecified)
			{
			// In the linear case, if the user did not specify the forces, then we
			// must generate the "f" array. Do this using linear interpolation of the
			// e array (which itself was generated using spline interpolation above).
			for (int i = 0; i < tablength; i++) {
			int im1 = i-1; if (im1 < 0) im1 += tablength;
			int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
			double dedx = (tb->e[ip1] - tb->e[im1]) / (2.0 * tb->delta);
			// (This is the average of the linear slopes on either side of the node.
			// Note that the nodes in the internal table are evenly spaced.)
			tb->f[i] = -dedx;
			}
			}

			// Fill the linear interpolation tables (de, df)
			for (int i = 0; i < tablength; i++) {
			int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
			tb->de[i] = tb->e[ip1] - tb->e[i];
			tb->df[i] = tb->f[ip1] - tb->f[i];
			}
			} // if (tabstyle == LINEAR)

			cyc_spline(tb->phi, tb->e, tablength, TWOPI, tb->e2);
			if (! tb->f_unspecified)
			cyc_spline(tb->phi, tb->f, tablength, TWOPI, tb->f2);
			}


			/* ----------------------------------------------------------------------
			extract attributes from parameter line in table section
			format of line: N value NOF DEGREES RADIANS
			N is required, other params are optional
			------------------------------------------------------------------------- */

			void DihedralTable::param_extract(Table tb, char line)
			{
			//tb->theta0 = 180.0; <- equilibrium angles not supported
			tb->ninput = 0;
			tb->f_unspecified = false; //default
			tb->use_degrees = true; //default
			strcpy(checkU_fname, "");
			strcpy(checkF_fname, "");

			char *word = strtok(line," \t\n\r\f");
			while (word) {
			if (strcmp(word,"N") == 0) {
			word = strtok(NULL," \t\n\r\f");
			tb->ninput = atoi(word);
			}
			else if (strcmp(word,"NOF") == 0) {
			tb->f_unspecified = true;
			}
			else if ((strcmp(word,"DEGREES") == 0) \|\| (strcmp(word,"degrees") == 0)) {
			tb->use_degrees = true;
			}
			else if ((strcmp(word,"RADIANS") == 0) \|\| (strcmp(word,"radians") == 0)) {
			tb->use_degrees = false;
			}
			else if (strcmp(word,"CHECKU") == 0) {
			word = strtok(NULL," \t\n\r\f");
			strncpy(checkU_fname, word, MAXLINE);
			}
			else if (strcmp(word,"CHECKF") == 0) {
			word = strtok(NULL," \t\n\r\f");
			strncpy(checkF_fname, word, MAXLINE);
			}
			// COMMENTING OUT: equilibrium angles are not supported
			//else if (strcmp(word,"EQ") == 0) {
			// word = strtok(NULL," \t\n\r\f");
			// tb->theta0 = atof(word);
			//}
			else {
			string err_msg("Invalid keyword in dihedral angle table parameters");
			err_msg += string(" (") + string(word) + string(")");
			error->one(FLERR,err_msg.c_str());
			}
			word = strtok(NULL," \t\n\r\f");
			}

			if (tb->ninput == 0)
			error->one(FLERR,"Dihedral table parameters did not set N");

			} // DihedralTable::param_extract()


			/* ----------------------------------------------------------------------
			broadcast read-in table info from proc 0 to other procs
			this function communicates these values in Table:
			ninput,phifile,efile,ffile,
			f_unspecified,use_degrees
			------------------------------------------------------------------------- */

			void DihedralTable::bcast_table(Table *tb)
			{
			MPI_Bcast(&tb->ninput,1,MPI_INT,0,world);

			int me;
			MPI_Comm_rank(world,&me);
			if (me > 0) {
			memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
			memory->create(tb->efile,tb->ninput,"dihedral:efile");
			memory->create(tb->ffile,tb->ninput,"dihedral:ffile");
			}

			MPI_Bcast(tb->phifile,tb->ninput,MPI_DOUBLE,0,world);
			MPI_Bcast(tb->efile,tb->ninput,MPI_DOUBLE,0,world);
			MPI_Bcast(tb->ffile,tb->ninput,MPI_DOUBLE,0,world);

			MPI_Bcast(&tb->f_unspecified,1,MPI_INT,0,world);
			MPI_Bcast(&tb->use_degrees,1,MPI_INT,0,world);

			// COMMENTING OUT: equilibrium angles are not supported
			//MPI_Bcast(&tb->theta0,1,MPI_DOUBLE,0,world);
			}


			namespace LAMMPS_NS {
			namespace DIHEDRAL_TABLE_NS {


			// -------------------------------------------------------------------
			// --------- The function was stolen verbatim from the ---------
			// --------- GNU Scientific Library (GSL, version 1.15) ---------
			// -------------------------------------------------------------------

			/* Author: Gerard Jungman */
			/* for description of method see [Engeln-Mullges + Uhlig, p. 96]
			*
			* diag[0] offdiag[0] 0 ..... offdiag[N-1]
			* offdiag[0] diag[1] offdiag[1] .....
			* 0 offdiag[1] diag[2]
			* 0 0 offdiag[2] .....
			* ... ...
			* offdiag[N-1] ...
			*
			*/
			// -- (A non-symmetric version of this function is also available.) --

			enum { //GSL status return codes.
			GSL_FAILURE = -1,
			GSL_SUCCESS = 0,
			GSL_ENOMEM = 8,
			GSL_EZERODIV = 12,
			GSL_EBADLEN = 19
			};


			int
			solve_cyc_tridiag(
			const double diag[], size_t d_stride,
			const double offdiag[], size_t o_stride,
			const double b[], size_t b_stride,
			double x[], size_t x_stride,
			size_t N)
			{
			int status = GSL_SUCCESS;
			double * delta = (double ) malloc (N sizeof (double));
			double * gamma = (double ) malloc (N sizeof (double));
			double * alpha = (double ) malloc (N sizeof (double));
			double * c = (double ) malloc (N sizeof (double));
			double * z = (double ) malloc (N sizeof (double));

			if (delta == 0 \|\| gamma == 0 \|\| alpha == 0 \|\| c == 0 \|\| z == 0) {
			cerr << "Internal Cyclic Spline Error: failed to allocate working space\n";
			exit(GSL_ENOMEM);
			}
			else
			{
			size_t i, j;
			double sum = 0.0;

			/* factor */

			if (N == 1)
			{
			x[0] = b[0] / diag[0];
			return GSL_SUCCESS;
			}

			alpha[0] = diag[0];
			gamma[0] = offdiag[0] / alpha[0];
			delta[0] = offdiag[o_stride * (N-1)] / alpha[0];

			if (alpha[0] == 0) {
			status = GSL_EZERODIV;
			}

			for (i = 1; i < N - 2; i++)
			{
			alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1];
			gamma[i] = offdiag[o_stride * i] / alpha[i];
			delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i];
			if (alpha[i] == 0) {
			status = GSL_EZERODIV;
			}
			}

			for (i = 0; i < N - 2; i++)
			{
			sum += alpha[i] * delta[i] * delta[i];
			}

			alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];

			gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];

			alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2];

			/* update */
			z[0] = b[0];
			for (i = 1; i < N - 1; i++)
			{
			z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1];
			}
			sum = 0.0;
			for (i = 0; i < N - 2; i++)
			{
			sum += delta[i] * z[i];
			}
			z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
			for (i = 0; i < N; i++)
			{
			c[i] = z[i] / alpha[i];
			}

			/* backsubstitution */
			x[x_stride * (N - 1)] = c[N - 1];
			x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)];
			if (N >= 3)
			{
			for (i = N - 3, j = 0; j <= N - 3; j++, i--)
			{
			x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)];
			}
			}
			}

			if (z != 0)
			free (z);
			if (c != 0)
			free (c);
			if (alpha != 0)
			free (alpha);
			if (gamma != 0)
			free (gamma);
			if (delta != 0)
			free (delta);

			if (status == GSL_EZERODIV) {
			cerr <<"Internal Cyclic Spline Error: Matrix must be positive definite.\n";
			exit(GSL_EZERODIV);
			}

			return status;
			} //solve_cyc_tridiag()



			/* ----------------------------------------------------------------------
			spline and splint routines modified from Numerical Recipes
			------------------------------------------------------------------------- */

			void cyc_spline(double const *xa,
			double const *ya,
			int n,
			double period,
			double *y2a)
			{

			double *diag = new double[n];
			double *offdiag = new double[n];
			double *rhs = new double[n];
			double xa_im1, xa_ip1;

			// In the cyclic case, there are n equations with n unknows.
			// The for loop sets up the equations we need to solve.
			// Later we invoke the GSL tridiagonal matrix solver to solve them.

			for(int i=0; i < n; i++) {

			// I have to lookup xa[i+1] and xa[i-1]. This gets tricky because of
			// periodic boundary conditions. We handle that now.
			int im1 = i-1;
			if (im1<0) {
			im1 += n;
			xa_im1 = xa[im1] - period;
			}
			else
			xa_im1 = xa[im1];

			int ip1 = i+1;
			if (ip1>=n) {
			ip1 -= n;
			xa_ip1 = xa[ip1] + period;
			}
			else
			xa_ip1 = xa[ip1];

			// Recall that we want to find the y2a[] parameters (there are n of them).
			// To solve for them, we have a linear equation with n unknowns
			// (in the cyclic case that is). For details, the non-cyclic case is
			// explained in equation 3.3.7 in Numerical Recipes in C, p. 115.
			diag[i] = (xa_ip1 - xa_im1) / 3.0;
			offdiag[i] = (xa_ip1 - xa[i]) / 6.0;
			rhs[i] = ((ya[ip1] - ya[i]) / (xa_ip1 - xa[i])) -
			((ya[i] - ya[im1]) / (xa[i] - xa_im1));
			}

			// Ordinarily we would have to invert a large matrix (with potentially
			// thousands of rows and columns). However because this matix happens
			// to be tridiagonal (and cyclic), we can use the following cheap method:
			solve_cyc_tridiag(diag, 1,
			offdiag, 1,
			rhs, 1,
			y2a, 1,
			n);

			delete [] diag;
			delete [] offdiag;
			delete [] rhs;

			} // cyc_spline()



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

			// cyc_splint(): Evaluates a spline at position x, with n control
			// points located at xa[], ya[], and with parameters y2a[]
			// The xa[] must be monotonically increasing and their
			// range should not exceed period (ie xa[n-1] < xa[0] + period).
			// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)]
			// "period" is typically 2*PI.
			double cyc_splint(double const *xa,
			double const *ya,
			double const *y2a,
			int n,
			double period,
			double x)
			{
			int klo = -1;
			int khi = n;
			int k;
			double xlo = xa[n-1] - period;
			double xhi = xa[0] + period;
			assert((xlo <= x) && (x <= xhi));
			while (khi-klo > 1) {
			k = (khi+klo) >> 1; //(k=(khi+klo)/2)
			if (xa[k] > x) {
			khi = k;
			xhi = xa[k];
			}
			else {
			klo = k;
			xlo = xa[k];
			}
			}
			assert((xlo <= x) && (x <= xhi));

			if (khi == n) khi = 0;
			if (klo ==-1) klo = n-1;

			double h = xhi-xlo;
			double a = (xhi-x) / h;
			double b = (x-xlo) / h;
			double y = aya[klo] + bya[khi] +
			((aaa-a)y2a[klo] + (bbb-b)y2a[khi]) * (h*h)/6.0;

			return y;

			} // cyc_splint()



			// cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
			// with n control points at xa[], ya[], with parameters y2a[].
			// The xa[] must be monotonically increasing and their
			// range should not exceed period (ie xa[n-1] < xa[0] + period).
			// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)]
			// "period" is typically 2*PI.
			double cyc_splintD(double const *xa,
			double const *ya,
			double const *y2a,
			int n,
			double period,
			double x)
			{
			int klo = -1;
			int khi = n; // (not n-1)
			int k;
			double xlo = xa[n-1] - period;
			double xhi = xa[0] + period;
			assert((xlo <= x) && (x <= xhi));
			while (khi-klo > 1) {
			k = (khi+klo) >> 1; //(k=(khi+klo)/2)
			if (xa[k] > x) {
			khi = k;
			xhi = xa[k];
			}
			else {
			klo = k;
			xlo = xa[k];
			}
			}
			assert((xlo <= x) && (x <= xhi));

			if (khi == n) khi = 0;
			if (klo ==-1) klo = n-1;

			double yhi = ya[khi];
			double ylo = ya[klo];
			double h = xhi-xlo;
			double g = yhi-ylo;
			double a = (xhi-x) / h;
			double b = (x-xlo) / h;
			// Formula below taken from equation 3.3.5 of "numerical recipes in c"
			// "yD" = the derivative of y
			double yD = g/h - ( (3.0aa-1.0)y2a[klo] - (3.0bb-1.0)y2a[khi] ) * h/6.0;
			// For rerefence: y = aylo + byhi +
			// ((aaa-a)y2a[klo] + (bbb-b)y2a[khi]) * (h*h)/6.0;

			return yD;

			} // cyc_splintD()



			} //namespace DIHEDRAL_TABLE_NS
			} //namespace LAMMPS_NS

lammps/src/USER-MISC/dihedral_table.cpp56933a0bd14alammm-master-deprecated

dihedral_table.cppView Options

lammps/src/USER-MISC/dihedral_table.cpp
56933a0bd14alammm-master-deprecated

dihedral_table.cpp
View Options