diff --git a/src/USER-MISC/Install.sh b/src/USER-MISC/Install.sh index d8e5b17d9f..c038252d60 100644 --- a/src/USER-MISC/Install.sh +++ b/src/USER-MISC/Install.sh @@ -10,12 +10,19 @@ if (test $1 = 1) then cp angle_cosine_shift.cpp .. cp angle_cosine_shift_exp.cpp .. cp angle_dipole.cpp .. + cp angle_fourier.cpp .. + cp angle_fourier_simple.cpp .. + cp angle_quartic.cpp .. cp bond_harmonic_shift.cpp .. cp bond_harmonic_shift_cut.cpp .. cp compute_ackland_atom.cpp .. cp compute_temp_rotate.cpp .. cp dihedral_cosine_shift_exp.cpp .. + cp dihedral_fourier.cpp .. + cp dihedral_nharmonic.cpp .. + cp dihedral_quadratic.cpp .. cp dihedral_table.cpp .. + cp improper_fourier.cpp .. cp fix_addtorque.cpp .. cp fix_imd.cpp .. cp fix_smd.cpp .. @@ -32,12 +39,19 @@ if (test $1 = 1) then cp angle_cosine_shift.h .. cp angle_cosine_shift_exp.h .. cp angle_dipole.h .. + cp angle_fourier.h .. + cp angle_fourier_simple.h .. + cp angle_quartic.h .. cp bond_harmonic_shift.h .. cp bond_harmonic_shift_cut.h .. cp compute_ackland_atom.h .. cp compute_temp_rotate.h .. cp dihedral_cosine_shift_exp.h .. + cp dihedral_fourier.h .. + cp dihedral_nharmonic.h .. + cp dihedral_quadratic.h .. cp dihedral_table.h .. + cp improper_fourier.h .. cp fix_addtorque.h .. cp fix_imd.h .. cp fix_smd.h .. @@ -55,13 +69,20 @@ elif (test $1 = 0) then rm -f ../angle_cosine_shift.cpp rm -f ../angle_cosine_shift_exp.cpp + rm -f ../angle_fourier.cpp + rm -f ../angle_fourier_simple.cpp rm -f ../angle_dipole.cpp + rm -f ../angle_quartic.cpp rm -f ../bond_harmonic_shift.cpp rm -f ../bond_harmonic_shift_cut.cpp rm -f ../compute_ackland_atom.cpp rm -f ../compute_temp_rotate.cpp rm -f ../dihedral_cosine_shift_exp.cpp + rm -f ../dihedral_fourier.cpp + rm -f ../dihedral_nharmonic.cpp + rm -f ../dihedral_quadratic.cpp rm -f ../dihedral_table.cpp + rm -f ../improper_fourier.cpp rm -f ../fix_addtorque.cpp rm -f ../fix_imd.cpp rm -f ../fix_smd.cpp @@ -79,12 +100,19 @@ elif (test $1 = 0) then rm -f ../angle_cosine_shift.h rm -f ../angle_cosine_shift_exp.h rm -f ../angle_dipole.h + rm -f ../angle_fourier.h + rm -f ../angle_fourier_simple.h + rm -f ../angle_quartic.h rm -f ../bond_harmonic_shift.h rm -f ../bond_harmonic_shift_cut.h rm -f ../compute_ackland_atom.h rm -f ../compute_temp_rotate.h rm -f ../dihedral_cosine_shift_exp.h + rm -f ../dihedral_fourier.h + rm -f ../dihedral_nharmonic.h + rm -f ../dihedral_quadratic.h rm -f ../dihedral_table.h + rm -f ../improper_fourier.h rm -f ../fix_addtorque.h rm -f ../fix_imd.h rm -f ../fix_smd.h diff --git a/src/USER-MISC/angle_fourier.cpp b/src/USER-MISC/angle_fourier.cpp new file mode 100644 index 0000000000..3f62b059ea --- /dev/null +++ b/src/USER-MISC/angle_fourier.cpp @@ -0,0 +1,268 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on angle_cosine_squared.cpp Naveen Michaud-Agrawal (Johns Hopkins U)] +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "angle_fourier.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleFourier::AngleFourier(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleFourier::~AngleFourier() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(C0); + memory->destroy(C1); + memory->destroy(C2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleFourier::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double term; + double rsq1,rsq2,r1,r2,c,c2,a,a11,a12,a22; + + eangle = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + + c2 = 2.0*c*c-1.0; + term = k[type]*(C0[type]+C1[type]*c+C2[type]*c2); + + if (eflag) eangle = term; + + a = k[type]*(C1[type]+4.0*C2[type]*c); + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + f1[0] = a11*delx1 + a12*delx2; + f1[1] = a11*dely1 + a12*dely2; + f1[2] = a11*delz1 + a12*delz2; + f3[0] = a22*delx2 + a12*delx1; + f3[1] = a22*dely2 + a12*dely1; + f3[2] = a22*delz2 + a12*delz1; + + // apply force to each of 3 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] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleFourier::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(k,n+1,"angle:k"); + memory->create(C0,n+1,"angle:C0"); + memory->create(C1,n+1,"angle:C1"); + memory->create(C2,n+1,"angle:C2"); + + memory->create(setflag,n+1,"angle:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void AngleFourier::coeff(int narg, char **arg) +{ + if (narg != 5) error->all(FLERR,"Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nangletypes,ilo,ihi); + + double k_one = force->numeric(arg[1]); + double C0_one = force->numeric(arg[2]); + double C1_one = force->numeric(arg[3]); + double C2_one = force->numeric(arg[4]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + C0[i] = C0_one; + C1[i] = C1_one; + C2[i] = C2_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleFourier::equilibrium_angle(int i) +{ + double ret=MY_PI; + if ( C2[i] != 0.0 ) { + ret = (C1[i]/4.0/C2[i]); + if ( abs(ret) <= 1.0 ) ret = acos(-ret); + } + return ret; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleFourier::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nangletypes,fp); + fwrite(&C0[1],sizeof(double),atom->nangletypes,fp); + fwrite(&C1[1],sizeof(double),atom->nangletypes,fp); + fwrite(&C2[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleFourier::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->nangletypes,fp); + fread(&C0[1],sizeof(double),atom->nangletypes,fp); + fread(&C1[1],sizeof(double),atom->nangletypes,fp); + fread(&C2[1],sizeof(double),atom->nangletypes,fp); + } + MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C0[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C1[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C2[1],atom->nangletypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- */ + +double AngleFourier::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double c2 = 2.0*c*c-1.0; + + double eng = k[type]*(C0[type]+C1[type]*c+C2[type]*c2); + return eng; +} diff --git a/src/USER-MISC/angle_fourier.h b/src/USER-MISC/angle_fourier.h new file mode 100644 index 0000000000..23ce20b4f4 --- /dev/null +++ b/src/USER-MISC/angle_fourier.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS + +AngleStyle(fourier,AngleFourier) + +#else + +#ifndef ANGLE_FOURIER_H +#define ANGLE_FOURIER_H + +#include "stdio.h" +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleFourier : public Angle { + public: + AngleFourier(class LAMMPS *); + virtual ~AngleFourier(); + virtual void compute(int, int); + void coeff(int, char **); + double equilibrium_angle(int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual double single(int, int, int, int); + + protected: + double *k,*C0,*C1,*C2; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/angle_fourier_simple.cpp b/src/USER-MISC/angle_fourier_simple.cpp new file mode 100644 index 0000000000..a67c721df4 --- /dev/null +++ b/src/USER-MISC/angle_fourier_simple.cpp @@ -0,0 +1,274 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on angle_cosine_quared.cpp Naveen Michaud-Agrawal (Johns Hopkins U)] +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "angle_fourier_simple.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleFourierSimple::AngleFourierSimple(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleFourierSimple::~AngleFourierSimple() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(C); + memory->destroy(N); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleFourierSimple::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double term,sgn; + double rsq1,rsq2,r1,r2,c,cn,th,nth,a,a11,a12,a22; + + eangle = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + + th = acos(c); + nth = N[type]*acos(c); + cn = cos(nth); + term = k[type]*(1.0+C[type]*cn); + + if (eflag) eangle = term; + + // handle sin(n th)/sin(th) singulatiries + + if ( abs(c)-1.0 > 0.0001 ) + { + a = k[type]*C[type]*N[type]*sin(nth)/sin(th); + } else { + if ( c >= 0.0 ) { + term = 1.0 - c; + sgn = 1.0; + } else { + term = 1.0 + c; + sgn = ( fmodf((float)(N[type]),2.0) == 0.0f )?-1.0:1.0; + } + a = N[type]+N[type]*(1.0-N[type]*N[type])*term/3.0; + a = k[type]*C[type]*N[type]*(double)(sgn)*a; + } + + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + f1[0] = a11*delx1 + a12*delx2; + f1[1] = a11*dely1 + a12*dely2; + f1[2] = a11*delz1 + a12*delz2; + f3[0] = a22*delx2 + a12*delx1; + f3[1] = a22*dely2 + a12*dely1; + f3[2] = a22*delz2 + a12*delz1; + + // apply force to each of 3 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] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleFourierSimple::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(k,n+1,"angle:k"); + memory->create(C,n+1,"angle:C"); + memory->create(N,n+1,"angle:N"); + + memory->create(setflag,n+1,"angle:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void AngleFourierSimple::coeff(int narg, char **arg) +{ + if (narg != 4) error->all(FLERR,"Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nangletypes,ilo,ihi); + + double k_one = force->numeric(arg[1]); + double C_one = force->numeric(arg[2]); + double N_one = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + C[i] = C_one; + N[i] = N_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleFourierSimple::equilibrium_angle(int i) +{ + return (MY_PI/N[i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleFourierSimple::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nangletypes,fp); + fwrite(&C[1],sizeof(double),atom->nangletypes,fp); + fwrite(&N[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleFourierSimple::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->nangletypes,fp); + fread(&C[1],sizeof(double),atom->nangletypes,fp); + fread(&N[1],sizeof(double),atom->nangletypes,fp); + } + MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&N[1],atom->nangletypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- */ + +double AngleFourierSimple::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double cn = cos(N[type]*acos(c)); + + double eng = k[type]*(1.0+C[type]*cn); + return eng; +} diff --git a/src/USER-MISC/angle_fourier_simple.h b/src/USER-MISC/angle_fourier_simple.h new file mode 100644 index 0000000000..411f04fa13 --- /dev/null +++ b/src/USER-MISC/angle_fourier_simple.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS + +AngleStyle(fourier/simple,AngleFourierSimple) + +#else + +#ifndef ANGLE_FOURIER_SIMPLE_H +#define ANGLE_FOURIER_SIMPLE_H + +#include "stdio.h" +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleFourierSimple : public Angle { + public: + AngleFourierSimple(class LAMMPS *); + virtual ~AngleFourierSimple(); + virtual void compute(int, int); + void coeff(int, char **); + double equilibrium_angle(int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual double single(int, int, int, int); + + protected: + double *k,*C,*N; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/angle_quartic.cpp b/src/USER-MISC/angle_quartic.cpp new file mode 100644 index 0000000000..ba310de67d --- /dev/null +++ b/src/USER-MISC/angle_quartic.cpp @@ -0,0 +1,277 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on angle_harmonic.cpp] +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "angle_quartic.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleQuartic::AngleQuartic(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleQuartic::~AngleQuartic() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k2); + memory->destroy(k3); + memory->destroy(k4); + memory->destroy(theta0); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleQuartic::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double dtheta,dtheta2,dtheta3,dtheta4,tk; + double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22; + + eangle = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + // force & energy + + dtheta = acos(c) - theta0[type]; + dtheta2 = dtheta * dtheta; + dtheta3 = dtheta2 * dtheta; + tk = 2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3; + + if (eflag) { + dtheta4 = dtheta3 * dtheta; + eangle = k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4; + } + + a = -2.0 * tk * s; + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + f1[0] = a11*delx1 + a12*delx2; + f1[1] = a11*dely1 + a12*dely2; + f1[2] = a11*delz1 + a12*delz2; + f3[0] = a22*delx2 + a12*delx1; + f3[1] = a22*dely2 + a12*dely1; + f3[2] = a22*delz2 + a12*delz1; + + // apply force to each of 3 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] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleQuartic::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(k2,n+1,"angle:k2"); + memory->create(k3,n+1,"angle:k3"); + memory->create(k4,n+1,"angle:k4"); + memory->create(theta0,n+1,"angle:theta0"); + + memory->create(setflag,n+1,"angle:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void AngleQuartic::coeff(int narg, char **arg) +{ + if (narg != 5) error->all(FLERR,"Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nangletypes,ilo,ihi); + + double theta0_one = force->numeric(arg[1]); + double k2_one = force->numeric(arg[2]); + double k3_one = force->numeric(arg[3]); + double k4_one = force->numeric(arg[4]); + + // convert theta0 from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k2[i] = k2_one; + k3[i] = k3_one; + k4[i] = k4_one; + theta0[i] = theta0_one/180.0 * MY_PI; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleQuartic::equilibrium_angle(int i) +{ + return theta0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleQuartic::write_restart(FILE *fp) +{ + fwrite(&k2[1],sizeof(double),atom->nangletypes,fp); + fwrite(&k3[1],sizeof(double),atom->nangletypes,fp); + fwrite(&k4[1],sizeof(double),atom->nangletypes,fp); + fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleQuartic::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k2[1],sizeof(double),atom->nangletypes,fp); + fread(&k3[1],sizeof(double),atom->nangletypes,fp); + fread(&k4[1],sizeof(double),atom->nangletypes,fp); + fread(&theta0[1],sizeof(double),atom->nangletypes,fp); + } + MPI_Bcast(&k2[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&k3[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&k4[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- */ + +double AngleQuartic::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double dtheta = acos(c) - theta0[type]; + double dtheta2 = dtheta * dtheta; + double dtheta3 = dtheta2 * dtheta; + double dtheta4 = dtheta3 * dtheta; + double tk = 2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3; + return k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4; +} diff --git a/src/USER-MISC/angle_quartic.h b/src/USER-MISC/angle_quartic.h new file mode 100644 index 0000000000..751c0a7abe --- /dev/null +++ b/src/USER-MISC/angle_quartic.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS + +AngleStyle(quartic,AngleQuartic) + +#else + +#ifndef LMP_ANGLE_QUARTIC_H +#define LMP_ANGLE_QUARTIC_H + +#include "stdio.h" +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleQuartic : public Angle { + public: + AngleQuartic(class LAMMPS *); + virtual ~AngleQuartic(); + virtual void compute(int, int); + void coeff(int, char **); + double equilibrium_angle(int); + void write_restart(FILE *); + void read_restart(FILE *); + double single(int, int, int, int); + + protected: + double *k2, *k3, *k4, *theta0; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/dihedral_fourier.cpp b/src/USER-MISC/dihedral_fourier.cpp new file mode 100644 index 0000000000..43438d42f3 --- /dev/null +++ b/src/USER-MISC/dihedral_fourier.cpp @@ -0,0 +1,411 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on dihedral_charmm.cpp Paul Crozier (SNL) ] +------------------------------------------------------------------------- */ + +#include "lmptype.h" +#include "mpi.h" +#include "math.h" +#include "stdlib.h" +#include "dihedral_fourier.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "domain.h" +#include "force.h" +#include "pair.h" +#include "update.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define TOLERANCE 0.05 + +/* ---------------------------------------------------------------------- */ + +DihedralFourier::DihedralFourier(LAMMPS *lmp) : Dihedral(lmp) {} + +/* ---------------------------------------------------------------------- */ + +DihedralFourier::~DihedralFourier() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(nterms); + + for (int i=1; i<= atom->ndihedraltypes; i++) { + if ( k[i] ) delete [] k[i]; + if ( multiplicity[i] ) delete [] multiplicity[i]; + if ( shift[i] ) delete [] shift[i]; + if ( cos_shift[i] ) delete [] cos_shift[i]; + if ( sin_shift[i] ) delete [] sin_shift[i]; + } + delete [] k; + delete [] multiplicity; + delete [] shift; + delete [] cos_shift; + delete [] sin_shift; + + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralFourier::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,i,j,m,n,type; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double edihedral,f1[3],f2[3],f3[3],f4[3]; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; + double df,df1_,ddf1_,fg,hg,fga,hgb,gaa,gbb; + double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; + double c,s,p,p_,sx2,sy2,sz2; + int itype,jtype; + double delx,dely,delz,rsq,r2inv,r6inv; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *atomtype = atom->type; + int **dihedrallist = neighbor->dihedrallist; + int ndihedrallist = neighbor->ndihedrallist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + double qqrd2e = force->qqrd2e; + + 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]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + domain->minimum_image(vb1x,vb1y,vb1z); + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + domain->minimum_image(vb2x,vb2y,vb2z); + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + domain->minimum_image(vb2xm,vb2ym,vb2zm); + + // 3rd bond + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + domain->minimum_image(vb3x,vb3y,vb3z); + + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; + + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rg = sqrt(rgsq); + + rginv = ra2inv = rb2inv = 0.0; + if (rg > 0) rginv = 1.0/rg; + if (rasq > 0) ra2inv = 1.0/rasq; + if (rbsq > 0) rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); + + c = (ax*bx + ay*by + az*bz)*rabinv; + s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) { + int me; + MPI_Comm_rank(world,&me); + if (screen) { + char str[128]; + sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d", + me,update->ntimestep, + atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]); + error->warning(FLERR,str,0); + fprintf(screen," 1st atom: %d %g %g %g\n", + me,x[i1][0],x[i1][1],x[i1][2]); + fprintf(screen," 2nd atom: %d %g %g %g\n", + me,x[i2][0],x[i2][1],x[i2][2]); + fprintf(screen," 3rd atom: %d %g %g %g\n", + me,x[i3][0],x[i3][1],x[i3][2]); + fprintf(screen," 4th atom: %d %g %g %g\n", + me,x[i4][0],x[i4][1],x[i4][2]); + } + } + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force and energy + // p = sum(i=1,nterms) k_i*(1+cos(n_i*phi-d_i) + // dp = dp / dphi + edihedral = 0.0; + df = 0.0; + for (j=0; jndihedraltypes; + + memory->create(nterms,n+1,"dihedral:nterms"); + k = new double * [n+1]; + multiplicity = new int * [n+1]; + shift = new int * [n+1]; + cos_shift = new double * [n+1]; + sin_shift = new double * [n+1]; + for (int i = 1; i <= n; i++) { + k[i] = cos_shift[i] = sin_shift[i] = 0; + multiplicity[i] = shift[i] = 0; + } + + memory->create(setflag,n+1,"dihedral:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void DihedralFourier::coeff(int narg, char **arg) +{ + if (narg < 4) error->all(FLERR,"Incorrect args for dihedral coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi); + + // require integer values of shift for backwards compatibility + // arbitrary phase angle shift could be allowed, but would break + // backwards compatibility and is probably not needed + + double k_one; + int multiplicity_one; + int shift_one; + int nterms_one = force->inumeric(arg[1]); + + if (nterms_one < 1) + error->all(FLERR,"Incorrect number of terms arg for dihedral coefficients"); + + if (2+3*nterms_one < narg) + error->all(FLERR,"Incorrect number of arguments for dihedral coefficients"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + nterms[i] = nterms_one; + k[i] = new double [nterms_one]; + multiplicity[i] = new int [nterms_one]; + shift[i] = new int [nterms_one]; + cos_shift[i] = new double [nterms_one]; + sin_shift[i] = new double [nterms_one]; + for (int j = 0; jnumeric(arg[offset+1]); + multiplicity_one = force->inumeric(arg[offset+2]); + shift_one = force->inumeric(arg[offset+3]); + k[i][j] = k_one; + multiplicity[i][j] = multiplicity_one; + shift[i][j] = shift_one; + cos_shift[i][j] = cos(MY_PI*shift_one/180.0); + sin_shift[i][j] = sin(MY_PI*shift_one/180.0); + } + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void DihedralFourier::write_restart(FILE *fp) +{ + + fwrite(&nterms[1],sizeof(int),atom->ndihedraltypes,fp); + for(int i = 1; i <= atom->ndihedraltypes; i++) { + fwrite(k[i],sizeof(double),nterms[i],fp); + fwrite(multiplicity[i],sizeof(int),nterms[i],fp); + fwrite(shift[i],sizeof(int),nterms[i],fp); + } + +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void DihedralFourier::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) + fread(&nterms[1],sizeof(int),atom->ndihedraltypes,fp); + + MPI_Bcast(&nterms[1],atom->ndihedraltypes,MPI_INT,0,world); + + // allocate + for (int i=1; i<=atom->ndihedraltypes; i++) { + k[i] = new double [nterms[i]]; + multiplicity[i] = new int [nterms[i]]; + shift[i] = new int [nterms[i]]; + cos_shift[i] = new double [nterms[i]]; + sin_shift[i] = new double [nterms[i]]; + } + + if (comm->me == 0) { + for (int i=1; i<=atom->ndihedraltypes; i++) { + fread(k[i],sizeof(double),nterms[i],fp); + fread(multiplicity[i],sizeof(int),nterms[i],fp); + fread(shift[i],sizeof(int),nterms[i],fp); + } + } + + for (int i=1; i<=atom->ndihedraltypes; i++) { + MPI_Bcast(k[i],nterms[i],MPI_DOUBLE,0,world); + MPI_Bcast(multiplicity[i],nterms[i],MPI_INT,0,world); + MPI_Bcast(shift[i],nterms[i],MPI_INT,0,world); + } + + for (int i=1; i <= atom->ndihedraltypes; i++) { + setflag[i] = 1; + for (int j = 0; j < nterms[i]; j++) { + cos_shift[i][j] = cos(MY_PI*shift[i][j]/180.0); + sin_shift[i][j] = sin(MY_PI*shift[i][j]/180.0); + } + } +} + diff --git a/src/USER-MISC/dihedral_fourier.h b/src/USER-MISC/dihedral_fourier.h new file mode 100644 index 0000000000..0c5dd863a3 --- /dev/null +++ b/src/USER-MISC/dihedral_fourier.h @@ -0,0 +1,49 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef DIHEDRAL_CLASS + +DihedralStyle(fourier,DihedralFourier) + +#else + +#ifndef LMP_DIHEDRAL_FOURIER_H +#define LMP_DIHEDRAL_FOURIER_H + +#include "stdio.h" +#include "dihedral.h" + +namespace LAMMPS_NS { + +class DihedralFourier : public Dihedral { + public: + DihedralFourier(class LAMMPS *); + virtual ~DihedralFourier(); + virtual void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + protected: + double **k,**cos_shift,**sin_shift; + int **multiplicity,**shift; + int *nterms; + int implicit,weightflag; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/dihedral_nharmonic.cpp b/src/USER-MISC/dihedral_nharmonic.cpp new file mode 100644 index 0000000000..1b0efb75d8 --- /dev/null +++ b/src/USER-MISC/dihedral_nharmonic.cpp @@ -0,0 +1,335 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on dihedral_multi_harmonic.cpp Mathias Puetz (SNL) and friends ] +------------------------------------------------------------------------- */ + +#include "lmptype.h" +#include "math.h" +#include "stdlib.h" +#include "dihedral_nharmonic.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define TOLERANCE 0.05 +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp) {} + +/* ---------------------------------------------------------------------- */ + +DihedralNHarmonic::~DihedralNHarmonic() +{ + if (allocated) { + memory->destroy(setflag); + for (int i = 1; i <= atom->ndihedraltypes; i++) + delete [] a[i]; + delete [] a; + delete [] nterms; + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralNHarmonic::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,n,type; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double edihedral,f1[3],f2[3],f3[3],f4[3]; + double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; + double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; + double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; + double a33,a12,a13,a23,sx2,sy2,sz2; + double s2,sin2; + + edihedral = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + 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; + + 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]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) { + int me; + MPI_Comm_rank(world,&me); + if (screen) { + char str[128]; + sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d", + me,update->ntimestep, + atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]); + error->warning(FLERR,str); + fprintf(screen," 1st atom: %d %g %g %g\n", + me,x[i1][0],x[i1][1],x[i1][2]); + fprintf(screen," 2nd atom: %d %g %g %g\n", + me,x[i2][0],x[i2][1],x[i2][2]); + fprintf(screen," 3rd atom: %d %g %g %g\n", + me,x[i3][0],x[i3][1],x[i3][2]); + fprintf(screen," 4th atom: %d %g %g %g\n", + me,x[i4][0],x[i4][1],x[i4][2]); + } + } + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + // p = sum (i=1,n) a_i * c**(i-1) + // pd = dp/dc + p = this->a[type][0]; + pd = this->a[type][1]; + for (int i = 1; i < nterms[type]-1; i++) { + p += c * this->a[type][i]; + pd += c * static_cast(i+1) * this->a[type][i+1]; + c *= c; + } + p += c * this->a[type][nterms[type]-1]; + + if (eflag) edihedral = p; + + a = pd; + c = c * a; + s12 = s12 * a; + a11 = c*sb1*s1; + a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); + a33 = c*sb3*s2; + a12 = -r12c1*(c1mag*c*s1 + c2mag*s12); + a13 = -rb1*rb3*s12; + a23 = r12c2*(c2mag*c*s2 + c1mag*s12); + + sx2 = a12*vb1x + a22*vb2x + a23*vb3x; + sy2 = a12*vb1y + a22*vb2y + a23*vb3y; + sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + + f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; + f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; + f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + + f2[0] = -sx2 - f1[0]; + f2[1] = -sy2 - f1[1]; + f2[2] = -sz2 - f1[2]; + + f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; + f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; + f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // 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, + vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralNHarmonic::allocate() +{ + allocated = 1; + int n = atom->ndihedraltypes; + + memory->create(nterms,n+1,"dihedral:nt"); + a = new double * [n+1]; + + memory->create(setflag,n+1,"dihedral:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void DihedralNHarmonic::coeff(int narg, char **arg) +{ + if (narg < 4 ) error->all(FLERR,"Incorrect args for dihedral coefficients"); + + int n = force->inumeric(arg[1]); + if (narg != n + 2 ) error->all(FLERR,"Incorrect args for dihedral coefficients"); + + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + a[i] = new double [n]; + nterms[i] = n; + for (int j = 0; j < n; j++ ) { + a[i][j] = force->numeric(arg[2+i]); + setflag[i] = 1; + } + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void DihedralNHarmonic::write_restart(FILE *fp) +{ + fwrite(&nterms[1],sizeof(int),atom->ndihedraltypes,fp); + for(int i = 1; i <= atom->ndihedraltypes; i++) + fwrite(a[i],sizeof(double),nterms[i],fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void DihedralNHarmonic::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) + fread(&nterms[1],sizeof(int),atom->ndihedraltypes,fp); + + MPI_Bcast(&nterms[1],atom->ndihedraltypes,MPI_INT,0,world); + + // allocate + for(int i = 1; i <= atom->ndihedraltypes; i++) + a[i] = new double [nterms[i]]; + + if (comm->me == 0) { + for(int i = 1; i <= atom->ndihedraltypes; i++) + fread(a[i],sizeof(double),nterms[i],fp); + } + + for (int i = 1; i <= atom->ndihedraltypes; i++ ) + MPI_Bcast(a[i],nterms[i],MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; +} + diff --git a/src/USER-MISC/dihedral_nharmonic.h b/src/USER-MISC/dihedral_nharmonic.h new file mode 100644 index 0000000000..31a66352d6 --- /dev/null +++ b/src/USER-MISC/dihedral_nharmonic.h @@ -0,0 +1,47 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef DIHEDRAL_CLASS + +DihedralStyle(n/harmonic,DihedralNHarmonic) + +#else + +#ifndef DIHEDRAL_NHARMONIC_H +#define DIHEDRAL_NHARMONIC_H + +#include "stdio.h" +#include "dihedral.h" + +namespace LAMMPS_NS { + +class DihedralNHarmonic : public Dihedral { + public: + DihedralNHarmonic(class LAMMPS *); + ~DihedralNHarmonic(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + protected: + int *nterms; + double **a; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/dihedral_quadratic.cpp b/src/USER-MISC/dihedral_quadratic.cpp new file mode 100644 index 0000000000..e617be70bd --- /dev/null +++ b/src/USER-MISC/dihedral_quadratic.cpp @@ -0,0 +1,338 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on dihedral_helix.cpp Paul Crozier (SNL) ] +------------------------------------------------------------------------- */ + +#include "lmptype.h" +#include "math.h" +#include "stdlib.h" +#include "dihedral_quadratic.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define TOLERANCE 0.05 +#define SMALL 0.001 +#define SMALLER 0.00001 + +/* ---------------------------------------------------------------------- */ + +DihedralQuadratic::DihedralQuadratic(LAMMPS *lmp) : Dihedral(lmp) {} + +/* ---------------------------------------------------------------------- */ + +DihedralQuadratic::~DihedralQuadratic() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(phi0); + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralQuadratic::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,i,m,n,type; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double edihedral,f1[3],f2[3],f3[3],f4[3]; + double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; + double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; + double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; + double a33,a12,a13,a23,sx2,sy2,sz2; + double s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + + edihedral = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + 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; + + 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]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + domain->minimum_image(vb1x,vb1y,vb1z); + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + domain->minimum_image(vb2x,vb2y,vb2z); + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + domain->minimum_image(vb2xm,vb2ym,vb2zm); + + // 3rd bond + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + domain->minimum_image(vb3x,vb3y,vb3z); + + // c0 calculation + + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + cx = vb1y*vb2z - vb1z*vb2y; + cy = vb1z*vb2x - vb1x*vb2z; + cz = vb1x*vb2y - vb1y*vb2x; + cmag = sqrt(cx*cx + cy*cy + cz*cz); + dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) { + int me; + MPI_Comm_rank(world,&me); + if (screen) { + char str[128]; + sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d", + me,update->ntimestep, + atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]); + error->warning(FLERR,str,0); + fprintf(screen," 1st atom: %d %g %g %g\n", + me,x[i1][0],x[i1][1],x[i1][2]); + fprintf(screen," 2nd atom: %d %g %g %g\n", + me,x[i2][0],x[i2][1],x[i2][2]); + fprintf(screen," 3rd atom: %d %g %g %g\n", + me,x[i3][0],x[i3][1],x[i3][2]); + fprintf(screen," 4th atom: %d %g %g %g\n", + me,x[i4][0],x[i4][1],x[i4][2]); + } + } + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // force & energy + // p = k ( phi- phi0)^2 + // pd = dp/dc + + phi = acos(c); + if (dx < 0.0) phi *= -1.0; + si = sin(phi); + + double dphi = phi-phi0[type]; + p = k[type]*dphi; + if (fabs(si) < SMALLER) { + pd = - 2.0 * k[type]; + } else { + pd = - 2.0 * p / si; + } + p = p * dphi; + + if (eflag) edihedral = p; + + a = pd; + c = c * a; + s12 = s12 * a; + a11 = c*sb1*s1; + a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); + a33 = c*sb3*s2; + a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); + a13 = -rb1*rb3*s12; + a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + + sx2 = a12*vb1x + a22*vb2x + a23*vb3x; + sy2 = a12*vb1y + a22*vb2y + a23*vb3y; + sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + + f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; + f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; + f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + + f2[0] = -sx2 - f1[0]; + f2[1] = -sy2 - f1[1]; + f2[2] = -sz2 - f1[2]; + + f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; + f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; + f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // 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, + vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + } +} + +/* ---------------------------------------------------------------------- */ + +void DihedralQuadratic::allocate() +{ + allocated = 1; + int n = atom->ndihedraltypes; + + memory->create(k,n+1,"dihedral:k"); + memory->create(phi0,n+1,"dihedral:phi0"); + + memory->create(setflag,n+1,"dihedral:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void DihedralQuadratic::coeff(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR,"Incorrect args for dihedral coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi); + + double k_one = force->numeric(arg[1]); + double phi0_one= force->numeric(arg[2]); + + // require k >= 0 + if (k_one < 0.0) + error->all(FLERR,"Incorrect coefficient arg for dihedral coefficients"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + phi0[i] = phi0_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void DihedralQuadratic::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); + fwrite(&phi0[1],sizeof(double),atom->ndihedraltypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void DihedralQuadratic::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->ndihedraltypes,fp); + fread(&phi0[1],sizeof(int),atom->ndihedraltypes,fp); + } + MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); + MPI_Bcast(&phi0[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; +} + diff --git a/src/USER-MISC/dihedral_quadratic.h b/src/USER-MISC/dihedral_quadratic.h new file mode 100644 index 0000000000..d1f5d2128b --- /dev/null +++ b/src/USER-MISC/dihedral_quadratic.h @@ -0,0 +1,47 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef DIHEDRAL_CLASS + +DihedralStyle(quadratic,DihedralQuadratic) + +#else + +#ifndef LMP_DIHEDRAL_QUADRATIC_H +#define LMP_DIHEDRAL_QUADRATIC_H + +#include "stdio.h" +#include "dihedral.h" + +namespace LAMMPS_NS { + +class DihedralQuadratic : public Dihedral { + public: + DihedralQuadratic(class LAMMPS *); + ~DihedralQuadratic(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + protected: + double *k,*phi0; + int *sign,*multiplicity; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/improper_fourier.cpp b/src/USER-MISC/improper_fourier.cpp new file mode 100644 index 0000000000..81d76c66e4 --- /dev/null +++ b/src/USER-MISC/improper_fourier.cpp @@ -0,0 +1,344 @@ +/* ---------------------------------------------------------------------- + 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: Loukas D. Peristeras (Scienomics SARL) + [ based on improper_umbrella.cpp Tod A Pascal (Caltech) ] +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "improper_fourier.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define TOLERANCE 0.05 +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +ImproperFourier::ImproperFourier(LAMMPS *lmp) : Improper(lmp) {} + +/* ---------------------------------------------------------------------- */ + +ImproperFourier::~ImproperFourier() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(C0); + memory->destroy(C1); + memory->destroy(C2); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperFourier::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,n,type; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + int **improperlist = neighbor->improperlist; + int nimproperlist = neighbor->nimproperlist; + + for (n = 0; n < nimproperlist; n++) { + i1 = improperlist[n][0]; + i2 = improperlist[n][1]; + i3 = improperlist[n][2]; + i4 = improperlist[n][3]; + type = improperlist[n][4]; + + // 1st bond + + vb1x = x[i2][0] - x[i1][0]; + vb1y = x[i2][1] - x[i1][1]; + vb1z = x[i2][2] - x[i1][2]; + domain->minimum_image(vb1x,vb1y,vb1z); + + // 2nd bond + + vb2x = x[i3][0] - x[i1][0]; + vb2y = x[i3][1] - x[i1][1]; + vb2z = x[i3][2] - x[i1][2]; + domain->minimum_image(vb2x,vb2y,vb2z); + + // 3rd bond + + vb3x = x[i4][0] - x[i1][0]; + vb3y = x[i4][1] - x[i1][1]; + vb3z = x[i4][2] - x[i1][2]; + domain->minimum_image(vb3x,vb3y,vb3z); + + addone(i1,i2,i3,i4, type,evflag,eflag, + vb1x, vb1y, vb1z, + vb2x, vb2y, vb2z, + vb3x, vb3y, vb3z); + if ( all[type] ) { + addone(i1,i4,i2,i3, type,evflag,eflag, + vb3x, vb3y, vb3z, + vb1x, vb1y, vb1z, + vb2x, vb2y, vb2z); + addone(i1,i3,i4,i2, type,evflag,eflag, + vb2x, vb2y, vb2z, + vb3x, vb3y, vb3z, + vb1x, vb1y, vb1z); + } + } +} + +void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int &i4, + const int &type,const int &evflag,const int &eflag, + const double &vb1x, const double &vb1y, const double &vb1z, + const double &vb2x, const double &vb2y, const double &vb2z, + const double &vb3x, const double &vb3y, const double &vb3z) +{ + int n; + double eimproper,f1[3],f2[3],f3[3],f4[3]; + double domega,c,c2,a,s,projhfg,dhax,dhay,dhaz,dahx,dahy,dahz,cotphi; + double ax,ay,az,ra2,rh2,ra,rh,rar,rhr,arx,ary,arz,hrx,hry,hrz; + + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + eimproper = 0.0; + + // c0 calculation + // A = vb1 X vb2 is perpendicular to IJK plane + + ax = vb1y*vb2z-vb1z*vb2y; + ay = vb1z*vb2x-vb1x*vb2z; + az = vb1x*vb2y-vb1y*vb2x; + ra2 = ax*ax+ay*ay+az*az; + rh2 = vb3x*vb3x+vb3y*vb3y+vb3z*vb3z; + ra = sqrt(ra2); + rh = sqrt(rh2); + if (ra < SMALL) ra = SMALL; + if (rh < SMALL) rh = SMALL; + + rar = 1/ra; + rhr = 1/rh; + arx = ax*rar; + ary = ay*rar; + arz = az*rar; + hrx = vb3x*rhr; + hry = vb3y*rhr; + hrz = vb3z*rhr; + + c = arx*hrx+ary*hry+arz*hrz; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) { + int me; + MPI_Comm_rank(world,&me); + if (screen) { + char str[128]; + sprintf(str, + "Improper problem: %d " BIGINT_FORMAT " %d %d %d %d", + me,update->ntimestep, + atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]); + error->warning(FLERR,str,0); + fprintf(screen," 1st atom: %d %g %g %g\n", me,x[i1][0],x[i1][1],x[i1][2]); + fprintf(screen," 2nd atom: %d %g %g %g\n", me,x[i2][0],x[i2][1],x[i2][2]); + fprintf(screen," 3rd atom: %d %g %g %g\n", me,x[i3][0],x[i3][1],x[i3][2]); + fprintf(screen," 4th atom: %d %g %g %g\n", me,x[i4][0],x[i4][1],x[i4][2]); + } + } + + if (c > 1.0) s = 1.0; + if (c < -1.0) s = -1.0; + + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + cotphi = c/s; + + projhfg = (vb3x*vb1x+vb3y*vb1y+vb3z*vb1z) / + sqrt(vb1x*vb1x+vb1y*vb1y+vb1z*vb1z); + projhfg += (vb3x*vb2x+vb3y*vb2y+vb3z*vb2z) / + sqrt(vb2x*vb2x+vb2y*vb2y+vb2z*vb2z); + if (projhfg > 0.0) { + s *= -1.0; + cotphi *= -1.0; + } + + // force and energy + // E = k ( C0 + C1 cos(w) + C2 cos(w) + + c2 = 2.0*s*s-1.0; + if (eflag) eimproper = k[type]*(C0[type]+C1[type]*s+C2[type]*c2); + + // dhax = diffrence between H and A in X direction, etc + + a = k[type]*(C1[type]+4.0*C2[type]*s)*cotphi; + dhax = hrx-c*arx; + dhay = hry-c*ary; + dhaz = hrz-c*arz; + + dahx = arx-c*hrx; + dahy = ary-c*hry; + dahz = arz-c*hrz; + + f2[0] = (dhay*vb1z - dhaz*vb1y)*rar; + f2[1] = (dhaz*vb1x - dhax*vb1z)*rar; + f2[2] = (dhax*vb1y - dhay*vb1x)*rar; + + f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar; + f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar; + f3[2] = (-dhax*vb2y + dhay*vb2x)*rar; + + f4[0] = dahx*rhr; + f4[1] = dahy*rhr; + f4[2] = dahz*rhr; + + f1[0] = -(f2[0] + f3[0] + f4[0]); + f1[1] = -(f2[1] + f3[1] + f4[1]); + f1[2] = -(f2[2] + f3[2] + f4[2]); + + // apply force to each of 4 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]*a; + f[i1][1] += f1[1]*a; + f[i1][2] += f1[2]*a; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += f3[0]*a; + f[i2][1] += f3[1]*a; + f[i2][2] += f3[2]*a; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f2[0]*a; + f[i3][1] += f2[1]*a; + f[i3][2] += f2[2]*a; + } + + if (newton_bond || i4 < nlocal) { + f[i4][0] += f4[0]*a; + f[i4][1] += f4[1]*a; + f[i4][2] += f4[2]*a; + } + + if (evflag) + ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4, + vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); +} + +/* ---------------------------------------------------------------------- */ + +void ImproperFourier::allocate() +{ + allocated = 1; + int n = atom->nimpropertypes; + + memory->create(k,n+1,"improper:k"); + memory->create(C0,n+1,"improper:C0"); + memory->create(C1,n+1,"improper:C1"); + memory->create(C2,n+1,"improper:C2"); + memory->create(all,n+1,"improper:C2"); + + memory->create(setflag,n+1,"improper:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void ImproperFourier::coeff(int narg, char **arg) +{ + + if ( narg != 5 && narg != 6 ) error->all(FLERR,"Incorrect args for improper coefficients"); + + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nimpropertypes,ilo,ihi); + + double k_one = force->numeric(arg[1]); + double C0_one = force->numeric(arg[2]); + double C1_one = force->numeric(arg[3]); + double C2_one = force->numeric(arg[4]); + int all_one = 1; + if ( narg == 6 ) all_one = force->inumeric(arg[5]); + + // convert w0 from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + C0[i] = C0_one; + C1[i] = C1_one; + C2[i] = C2_one; + all[i] = all_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void ImproperFourier::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&C0[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&C1[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&C2[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&all[1],sizeof(int),atom->nimpropertypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void ImproperFourier::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->nimpropertypes,fp); + fread(&C0[1],sizeof(double),atom->nimpropertypes,fp); + fread(&C1[1],sizeof(double),atom->nimpropertypes,fp); + fread(&C2[1],sizeof(double),atom->nimpropertypes,fp); + fread(&all[1],sizeof(int),atom->nimpropertypes,fp); + } + MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C0[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C1[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&C2[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&all[1],atom->nimpropertypes,MPI_INT,0,world); + + for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1; +} diff --git a/src/USER-MISC/improper_fourier.h b/src/USER-MISC/improper_fourier.h new file mode 100644 index 0000000000..cfd47de135 --- /dev/null +++ b/src/USER-MISC/improper_fourier.h @@ -0,0 +1,51 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef IMPROPER_CLASS + +ImproperStyle(fourier,ImproperFourier) + +#else + +#ifndef LMP_IMPROPER_FOURIER_H +#define LMP_IMPROPER_FOURIER_H + +#include "stdio.h" +#include "improper.h" + +namespace LAMMPS_NS { + +class ImproperFourier : public Improper { + public: + ImproperFourier(class LAMMPS *); + ~ImproperFourier(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + protected: + double *k, *C0, *C1, *C2; + int *all; + void addone(const int &i1,const int &i2,const int &i3,const int &i4, + const int &type,const int &evflag,const int &eflag, + const double &vb1x, const double &vb1y, const double &vb1z, + const double &vb2x, const double &vb2y, const double &vb2z, + const double &vb3x, const double &vb3y, const double &vb3z); + void allocate(); +}; + +} + +#endif +#endif