From 67bcfb00e101522c70787db035ca907fcf90fd10 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 8 Aug 2011 22:55:05 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6647 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/angle_cosine_shift.cpp | 263 +++++++++++++++++++ src/USER-MISC/angle_cosine_shift.h | 52 ++++ src/USER-MISC/angle_cosine_shift_exp.cpp | 306 +++++++++++++++++++++++ src/USER-MISC/angle_cosine_shift_exp.h | 50 ++++ 4 files changed, 671 insertions(+) create mode 100644 src/USER-MISC/angle_cosine_shift.cpp create mode 100644 src/USER-MISC/angle_cosine_shift.h create mode 100644 src/USER-MISC/angle_cosine_shift_exp.cpp create mode 100644 src/USER-MISC/angle_cosine_shift_exp.h diff --git a/src/USER-MISC/angle_cosine_shift.cpp b/src/USER-MISC/angle_cosine_shift.cpp new file mode 100644 index 0000000000..4e0d11d849 --- /dev/null +++ b/src/USER-MISC/angle_cosine_shift.cpp @@ -0,0 +1,263 @@ +/* ---------------------------------------------------------------------- + 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: Carsten Svaneborg, science@zqex.dk +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "angle_cosineshift.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleCosineShift::AngleCosineShift(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleCosineShift::~AngleCosineShift() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(kcost); + memory->destroy(ksint); + memory->destroy(theta); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleCosineShift::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 rsq1,rsq2,r1,r2,c,s,cps,kcos,ksin,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]; + domain->minimum_image(delx1,dely1,delz1); + + 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]; + domain->minimum_image(delx2,dely2,delz2); + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // c = cosine of angle + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // C= sine of angle + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // force & energy + if (eflag) eangle = -k[type]-kcos*c-ksin*s; + + kcos=kcost[type]; + ksin=ksint[type]; + cps = c/s; // NOTE absorbed one c + + a11 = (-kcos +ksin*cps )*c/ rsq1; + a12 = ( kcos -ksin*cps ) / (r1*r2); + a22 = (-kcos +ksin*cps )*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 AngleCosineShift::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(k ,n+1,"Angle:k"); + memory->create(ksint ,n+1,"Angle:ksint"); + memory->create(kcost ,n+1,"Angle:kcost"); + memory->create(theta ,n+1,"Angle:theta"); + memory->create(setflag,n+1, "Angle:setflag"); + + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void AngleCosineShift::coeff(int narg, char **arg) +{ + if (narg != 3) error->all("Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nangletypes,ilo,ihi); + + double umin = force->numeric(arg[1]); + double theta0 = force->numeric(arg[2]); + +// k=Umin/2 + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = umin/2; + kcost[i] = umin/2*cos(theta0*3.14159265/180); + ksint[i] = umin/2*sin(theta0*3.14159265/180); + theta[i] = theta0*3.14159265/180; + + setflag[i] = 1; + count++; + } + + if (count == 0) error->all("Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleCosineShift::equilibrium_angle(int i) +{ + return theta[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleCosineShift::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nangletypes,fp); + fwrite(&kcost[1],sizeof(double),atom->nangletypes,fp); + fwrite(&ksint[1],sizeof(double),atom->nangletypes,fp); + fwrite(&theta[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleCosineShift::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) + { + fread(&k[1],sizeof(double),atom->nangletypes,fp); + fread(&kcost[1],sizeof(double),atom->nangletypes,fp); + fread(&ksint[1],sizeof(double),atom->nangletypes,fp); + fread(&theta[1],sizeof(double),atom->nangletypes,fp); + } + MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&kcost[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&ksint[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&theta[1],atom->nangletypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- */ + +double AngleCosineShift::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 s=sqrt(1.0-c*c); + + return -k[type]-kcost[type]*c-ksint[type]*s; +} diff --git a/src/USER-MISC/angle_cosine_shift.h b/src/USER-MISC/angle_cosine_shift.h new file mode 100644 index 0000000000..ce44ea3045 --- /dev/null +++ b/src/USER-MISC/angle_cosine_shift.h @@ -0,0 +1,52 @@ +/* ---------------------------------------------------------------------- + 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(cosine/shift,AngleCosineShift) + +#else + +#ifndef LMP_ANGLE_COSINE_SHIFT_H +#define LMP_ANGLE_COSINE_SHIFT_H + +#include "stdio.h" +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleCosineShift : public Angle { + public: + AngleCosineShift(class LAMMPS *); + ~AngleCosineShift(); + 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); + + private: + double *k; + double *a; + double *theta; + double *ksint; + double *kcost; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/angle_cosine_shift_exp.cpp b/src/USER-MISC/angle_cosine_shift_exp.cpp new file mode 100644 index 0000000000..f4a95ca7b0 --- /dev/null +++ b/src/USER-MISC/angle_cosine_shift_exp.cpp @@ -0,0 +1,306 @@ +/* ---------------------------------------------------------------------- + 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: Carsten Svaneborg, science@zqex.dk +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "angle_cosineshiftexp.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleCosineShiftExp::AngleCosineShiftExp(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleCosineShiftExp::~AngleCosineShiftExp() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(umin); + memory->destroy(a); + memory->destroy(opt1); + memory->destroy(cost); + memory->destroy(sint); + memory->destroy(theta0); + memory->destroy(doExpansion); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleCosineShiftExp::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3],ff; + double rsq1,rsq2,r1,r2,c,s,cc,ss,a11,a12,a22; + double exp2,aa,uumin,cccpsss,cssmscc; + + 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]; + domain->minimum_image(delx1,dely1,delz1); + + 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]; + domain->minimum_image(delx2,dely2,delz2); + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // c = cosine of angle + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // C= sine of angle + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // force & energy + + aa=a[type]; + uumin=umin[type]; + + cccpsss = c*cost[type]+s*sint[type]; + cssmscc = c*sint[type]-s*cost[type]; + + if (doExpansion[type]) + { // |a|<0.01 so use expansions relative precision <1e-5 +// std::cout << "Using expansion\n"; + if (eflag) eangle = -0.125*(1+cccpsss)*(4+aa*(cccpsss-1))*uumin; + ff=0.25*uumin*cssmscc*(2+aa*cccpsss)/s; + } + else + { +// std::cout << "Not using expansion\n"; + exp2=exp(0.5*aa*(1+cccpsss)); + if (eflag) eangle = opt1[type]*(1-exp2); + ff=0.5*a[type]*opt1[type]*exp2*cssmscc/s; + } + + a11 = ff*c/ rsq1; + a12 = -ff / (r1*r2); + a22 = ff*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 AngleCosineShiftExp::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(doExpansion, n+1, "angle:doExpansion"); + memory->create(umin , n+1, "angle:umin"); + memory->create(a , n+1, "angle:a"); + memory->create(sint , n+1, "angle:sint"); + memory->create(cost , n+1, "angle:cost"); + memory->create(opt1 , n+1, "angle:opt1"); + 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 type +------------------------------------------------------------------------- */ + +void AngleCosineShiftExp::coeff(int narg, char **arg) +{ + if (narg != 4) error->all("Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nangletypes,ilo,ihi); + + double umin_ = force->numeric(arg[1]); + double theta0_ = force->numeric(arg[2]); + double a_ = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + doExpansion[i]=(fabs(a_)<0.001); + umin[i] = umin_; + a[i] = a_; + cost[i] = cos(theta0_*3.14159265/180); + sint[i] = sin(theta0_*3.14159265/180); + theta0[i]= theta0_*3.14159265/180; + + if (!doExpansion[i]) opt1[i]=umin_/(exp(a_)-1); + + setflag[i] = 1; + count++; + } + + if (count == 0) error->all("Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleCosineShiftExp::equilibrium_angle(int i) +{ + return theta0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleCosineShiftExp::write_restart(FILE *fp) +{ + fwrite(&umin[1],sizeof(double),atom->nangletypes,fp); + fwrite(&a[1],sizeof(double),atom->nangletypes,fp); + fwrite(&cost[1],sizeof(double),atom->nangletypes,fp); + fwrite(&sint[1],sizeof(double),atom->nangletypes,fp); + fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleCosineShiftExp::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) + { + fread(&umin[1],sizeof(double),atom->nangletypes,fp); + fread(&a[1],sizeof(double),atom->nangletypes,fp); + fread(&cost[1],sizeof(double),atom->nangletypes,fp); + fread(&sint[1],sizeof(double),atom->nangletypes,fp); + fread(&theta0[1],sizeof(double),atom->nangletypes,fp); + } + MPI_Bcast(&umin[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&a[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&cost[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&sint[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; + doExpansion[i]=(fabs(a[i])<0.01); + if (!doExpansion[i]) opt1[i]=umin[i]/(exp(a[i])-1); + } +} + +/* ---------------------------------------------------------------------- */ + +double AngleCosineShiftExp::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 s=sqrt(1.0-c*c); + + double cccpsss=c*cost[type]+s*sint[type]; + double cssmscc=c*sint[type]-s*cost[type]; + + if (doExpansion[type]) + { + return -0.125*(1+cccpsss)*(4+a[type]*(cccpsss-1))*umin[type]; + } + else + { + return opt1[type]*(1-exp(0.5*a[type]*(1+cccpsss))); + } +} diff --git a/src/USER-MISC/angle_cosine_shift_exp.h b/src/USER-MISC/angle_cosine_shift_exp.h new file mode 100644 index 0000000000..c51ece7c44 --- /dev/null +++ b/src/USER-MISC/angle_cosine_shift_exp.h @@ -0,0 +1,50 @@ +/* ---------------------------------------------------------------------- + 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(cosine/shift/exp,AngleCosineShiftExp) +#else + +#ifndef LMP_ANGLE_COSINE_SHIFT_EXP_H +#define LMP_ANGLE_COSINE_SHIFT_EXP_H + +#include "stdio.h" +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleCosineShiftExp : public Angle { + public: + AngleCosineShiftExp(class LAMMPS *); + ~AngleCosineShiftExp(); + 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); + + private: + bool *doExpansion; + double *umin,*a,*opt1; + double *theta0; + double *sint; + double *cost; + + void allocate(); +}; + +} + +#endif +#endif