API tools faq
paste
Advertisement
Guest User

pair_lj_charmmfsw2_coul_long

a guest
Dec 1st, 2015
81
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 42.72 KB | None | 0 0
raw download clone embed print report
  1. /* ----------------------------------------------------------------------
  2.    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
  3.    http://lammps.sandia.gov, Sandia National Laboratories
  4.    Steve Plimpton, sjplimp@sandia.gov
  5.  
  6.    Copyright (2003) Sandia Corporation.  Under the terms of Contract
  7.    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
  8.    certain rights in this software.  This software is distributed under
  9.    the GNU General Public License.
  10.  
  11.    See the README file in the top-level LAMMPS directory.
  12. ------------------------------------------------------------------------- */
  13.  
  14. /* ----------------------------------------------------------------------
  15.    Contributing author: Paul Crozier (SNL)
  16. ------------------------------------------------------------------------- */
  17.  
  18. #include "math.h"
  19. #include "stdio.h"
  20. #include "stdlib.h"
  21. #include "string.h"
  22. #include "pair_lj_charmm_coul_long.h"
  23. #include "atom.h"
  24. #include "comm.h"
  25. #include "force.h"
  26. #include "kspace.h"
  27. #include "update.h"
  28. #include "integrate.h"
  29. #include "respa.h"
  30. #include "neighbor.h"
  31. #include "neigh_list.h"
  32. #include "neigh_request.h"
  33. #include "memory.h"
  34. #include "error.h"
  35.  
  36. using namespace LAMMPS_NS;
  37.  
  38. #define EWALD_F   1.12837917
  39. #define EWALD_P   0.3275911
  40. #define A1        0.254829592
  41. #define A2       -0.284496736
  42. #define A3        1.421413741
  43. #define A4       -1.453152027
  44. #define A5        1.061405429
  45.  
  46. /* ---------------------------------------------------------------------- */
  47.  
  48. PairLJCharmmCoulLong::PairLJCharmmCoulLong(LAMMPS *lmp) : Pair(lmp)
  49. {
  50.   respa_enable = 1;
  51.   ewaldflag = pppmflag = 1;
  52.   ftable = NULL;
  53.   implicit = 0;
  54.   mix_flag = ARITHMETIC;
  55.   writedata = 1;
  56. }
  57.  
  58. /* ---------------------------------------------------------------------- */
  59.  
  60. PairLJCharmmCoulLong::~PairLJCharmmCoulLong()
  61. {
  62.   if (!copymode) {
  63.     if (allocated) {
  64.       memory->destroy(setflag);
  65.       memory->destroy(cutsq);
  66.  
  67.       memory->destroy(epsilon);
  68.       memory->destroy(sigma);
  69.       memory->destroy(eps14);
  70.       memory->destroy(sigma14);
  71.       memory->destroy(lj1);
  72.       memory->destroy(lj2);
  73.       memory->destroy(lj3);
  74.       memory->destroy(lj4);
  75.       memory->destroy(lj14_1);
  76.       memory->destroy(lj14_2);
  77.       memory->destroy(lj14_3);
  78.       memory->destroy(lj14_4);
  79.       memory->destroy(offset);
  80.     }
  81.     if (ftable) free_tables();
  82.   }
  83. }
  84.  
  85. /* ---------------------------------------------------------------------- */
  86.  
  87. void PairLJCharmmCoulLong::compute(int eflag, int vflag)
  88. {
  89.   int i,j,ii,jj,inum,jnum,itype,jtype,itable;
  90.   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
  91.   double fraction,table;
  92.   double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
  93.   double grij,expm2,prefactor,t,erfc;
  94.   double philj,switch1,switch2;
  95.   int *ilist,*jlist,*numneigh,**firstneigh;
  96.   double rsq;
  97.  
  98.   evdwl = ecoul = 0.0;
  99.   if (eflag || vflag) ev_setup(eflag,vflag);
  100.   else evflag = vflag_fdotr = 0;
  101.  
  102.   double **x = atom->x;
  103.   double **f = atom->f;
  104.   double *q = atom->q;
  105.   int *type = atom->type;
  106.   int nlocal = atom->nlocal;
  107.   double *special_coul = force->special_coul;
  108.   double *special_lj = force->special_lj;
  109.   int newton_pair = force->newton_pair;
  110.   double qqrd2e = force->qqrd2e;
  111.  
  112.   inum = list->inum;
  113.   ilist = list->ilist;
  114.   numneigh = list->numneigh;
  115.   firstneigh = list->firstneigh;
  116.  
  117.   double cut_ljsq_6;
  118.   double cut_ljsq_3;
  119.   double cut_lj_innersq_6;
  120.   double cut_lj_innersq_3;
  121.   double k12;
  122.   double k6;
  123.   double r3inv;
  124.  
  125.   // only need to compute these once
  126.   cut_ljsq_6 = cut_ljsq*cut_ljsq*cut_ljsq;
  127.   cut_ljsq_3 = sqrt(cut_ljsq_6);
  128.   cut_lj_innersq_6 = cut_lj_innersq*cut_lj_innersq*cut_lj_innersq;
  129.   cut_lj_innersq_3 = sqrt(cut_lj_innersq_6);
  130.  
  131.  
  132.   // loop over neighbors of my atoms
  133.  
  134.   for (ii = 0; ii < inum; ii++) {
  135.     i = ilist[ii];
  136.     qtmp = q[i];
  137.     xtmp = x[i][0];
  138.     ytmp = x[i][1];
  139.     ztmp = x[i][2];
  140.     itype = type[i];
  141.     jlist = firstneigh[i];
  142.     jnum = numneigh[i];
  143.  
  144.     for (jj = 0; jj < jnum; jj++) {
  145.       j = jlist[jj];
  146.       factor_lj = special_lj[sbmask(j)];
  147.       factor_coul = special_coul[sbmask(j)];
  148.       j &= NEIGHMASK;
  149.  
  150.       delx = xtmp - x[j][0];
  151.       dely = ytmp - x[j][1];
  152.       delz = ztmp - x[j][2];
  153.       rsq = delx*delx + dely*dely + delz*delz;
  154.  
  155.       if (rsq < cut_bothsq) {
  156.         r2inv = 1.0/rsq;
  157.  
  158.         if (rsq < cut_coulsq) {
  159.           if (!ncoultablebits || rsq <= tabinnersq) {
  160.             r = sqrt(rsq);
  161.             grij = g_ewald * r;
  162.             expm2 = exp(-grij*grij);
  163.             t = 1.0 / (1.0 + EWALD_P*grij);
  164.             erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
  165.             prefactor = qqrd2e * qtmp*q[j]/r;
  166.             forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
  167.             if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
  168.           } else {
  169.             union_int_float_t rsq_lookup;
  170.             rsq_lookup.f = rsq;
  171.             itable = rsq_lookup.i & ncoulmask;
  172.             itable >>= ncoulshiftbits;
  173.             fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
  174.             table = ftable[itable] + fraction*dftable[itable];
  175.             forcecoul = qtmp*q[j] * table;
  176.             if (factor_coul < 1.0) {
  177.               table = ctable[itable] + fraction*dctable[itable];
  178.               prefactor = qtmp*q[j] * table;
  179.               forcecoul -= (1.0-factor_coul)*prefactor;
  180.             }
  181.           }
  182.         } else forcecoul = 0.0;
  183.  
  184. /* EVL modify
  185.         if (rsq < cut_ljsq) {
  186.           r6inv = r2inv*r2inv*r2inv;
  187.           jtype = type[j];
  188.           forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  189.           if (rsq > cut_lj_innersq) {
  190.             switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  191.               (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  192.             switch2 = 12.0*rsq * (cut_ljsq-rsq) *
  193.               (rsq-cut_lj_innersq) / denom_lj;
  194.             philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
  195.             forcelj = forcelj*switch1 + philj*switch2;
  196.           }
  197.         } else forcelj = 0.0;
  198. */
  199.         /* FORCE CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  200.         if (rsq <= cut_ljsq) {
  201.             r6inv = r2inv*r2inv*r2inv;
  202.             jtype = type[j];
  203.  
  204.           if (rsq <= cut_lj_innersq) {
  205.             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  206.  
  207.           } else {
  208.             r3inv = sqrt(r6inv);
  209.  
  210.             /* Steinbach and Brooks (1993) Eqs. 10b and 13 */
  211.             k12 =  lj3[itype][jtype] * cut_ljsq_6 / (cut_ljsq_6-cut_lj_innersq_6);
  212.             k6  = -lj4[itype][jtype] * cut_ljsq_3 / (cut_ljsq_3-cut_lj_innersq_3);
  213.  
  214.             /* Neg. Gradient of Steinbach and Brooks Eq. 10b  */
  215.             forcelj = 0;
  216.             forcelj += 2*6* k12*r6inv*(r6inv-1/cut_ljsq_6);
  217.             forcelj += 2*3* k6 *r3inv*(r3inv-1/cut_ljsq_3);
  218.           }
  219.  
  220.         } else forcelj = 0.0;
  221.  
  222.  
  223.         fpair = (forcecoul + factor_lj*forcelj) * r2inv;
  224.  
  225.         f[i][0] += delx*fpair;
  226.         f[i][1] += dely*fpair;
  227.         f[i][2] += delz*fpair;
  228.         if (newton_pair || j < nlocal) {
  229.           f[j][0] -= delx*fpair;
  230.           f[j][1] -= dely*fpair;
  231.           f[j][2] -= delz*fpair;
  232.         }
  233.  
  234.         if (eflag) {
  235.           if (rsq < cut_coulsq) {
  236.             if (!ncoultablebits || rsq <= tabinnersq)
  237.               ecoul = prefactor*erfc;
  238.             else {
  239.               table = etable[itable] + fraction*detable[itable];
  240.               ecoul = qtmp*q[j] * table;
  241.             }
  242.             if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
  243.           } else ecoul = 0.0;
  244.  
  245. /* EVL modify
  246.           if (rsq < cut_ljsq) {
  247.             evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
  248.             if (rsq > cut_lj_innersq) {
  249.               switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  250.                 (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  251.               evdwl *= switch1;
  252.             }
  253.             evdwl *= factor_lj;
  254.           } else evdwl = 0.0;
  255. */
  256.  
  257.         /* ENERGY CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  258.           if (rsq <= cut_ljsq) {
  259.             if (rsq <= cut_lj_innersq) {
  260.               /* Steinbach and Brooks (1993) Eqs. 10a and 12 */
  261.               philj  = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
  262.               philj += (-lj3[itype][jtype] / (cut_ljsq_6*cut_lj_innersq_6));
  263.               philj += ( lj4[itype][jtype] / (cut_ljsq_3*cut_lj_innersq_3));
  264.             } else {
  265.               /* Steinbach and Brooks (1993) Eqs. 10b */
  266.               philj = 0;
  267.               philj += k12 * (r6inv - 1/cut_ljsq_6) * (r6inv - 1/cut_ljsq_6);
  268.               philj += k6  * (r3inv - 1/cut_ljsq_3) * (r3inv - 1/cut_ljsq_3);
  269.             }
  270.             evdwl = factor_lj*philj;
  271.           } else evdwl = 0.0;
  272.  
  273.  
  274.         }
  275.  
  276.         if (evflag) ev_tally(i,j,nlocal,newton_pair,
  277.                              evdwl,ecoul,fpair,delx,dely,delz);
  278.       }
  279.     }
  280.   }
  281.  
  282.   if (vflag_fdotr) virial_fdotr_compute();
  283. }
  284.  
  285. /* ---------------------------------------------------------------------- */
  286.  
  287. void PairLJCharmmCoulLong::compute_inner()
  288. {
  289.   int i,j,ii,jj,inum,jnum,itype,jtype;
  290.   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
  291.   double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
  292.   double rsw;
  293.   int *ilist,*jlist,*numneigh,**firstneigh;
  294.  
  295.   double **x = atom->x;
  296.   double **f = atom->f;
  297.   double *q = atom->q;
  298.   int *type = atom->type;
  299.   int nlocal = atom->nlocal;
  300.   double *special_coul = force->special_coul;
  301.   double *special_lj = force->special_lj;
  302.   int newton_pair = force->newton_pair;
  303.   double qqrd2e = force->qqrd2e;
  304.  
  305.   inum = listinner->inum;
  306.   ilist = listinner->ilist;
  307.   numneigh = listinner->numneigh;
  308.   firstneigh = listinner->firstneigh;
  309.  
  310.   double cut_out_on = cut_respa[0];
  311.   double cut_out_off = cut_respa[1];
  312.  
  313.   double cut_out_diff = cut_out_off - cut_out_on;
  314.   double cut_out_on_sq = cut_out_on*cut_out_on;
  315.   double cut_out_off_sq = cut_out_off*cut_out_off;
  316.  
  317.   double cut_ljsq_6;
  318.   double cut_ljsq_3;
  319.   double cut_lj_innersq_6;
  320.   double cut_lj_innersq_3;
  321.   double k12;
  322.   double k6;
  323.   double r3inv;
  324.  
  325.   // only need to compute these once
  326.   cut_ljsq_6 = cut_ljsq*cut_ljsq*cut_ljsq;
  327.   cut_ljsq_3 = sqrt(cut_ljsq_6);
  328.   cut_lj_innersq_6 = cut_lj_innersq*cut_lj_innersq*cut_lj_innersq;
  329.   cut_lj_innersq_3 = sqrt(cut_lj_innersq_6);
  330.  
  331.  
  332.  
  333.   // loop over neighbors of my atoms
  334.  
  335.   for (ii = 0; ii < inum; ii++) {
  336.     i = ilist[ii];
  337.     qtmp = q[i];
  338.     xtmp = x[i][0];
  339.     ytmp = x[i][1];
  340.     ztmp = x[i][2];
  341.     itype = type[i];
  342.     jlist = firstneigh[i];
  343.     jnum = numneigh[i];
  344.  
  345.     for (jj = 0; jj < jnum; jj++) {
  346.       j = jlist[jj];
  347.       factor_lj = special_lj[sbmask(j)];
  348.       factor_coul = special_coul[sbmask(j)];
  349.       j &= NEIGHMASK;
  350.  
  351.       delx = xtmp - x[j][0];
  352.       dely = ytmp - x[j][1];
  353.       delz = ztmp - x[j][2];
  354.       rsq = delx*delx + dely*dely + delz*delz;
  355.  
  356.       if (rsq < cut_out_off_sq) {
  357.         r2inv = 1.0/rsq;
  358.         forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
  359.         if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
  360.  
  361.  
  362.         /* FORCE CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  363.         if (rsq <= cut_ljsq) {
  364.             r6inv = r2inv*r2inv*r2inv;
  365.             jtype = type[j];
  366.  
  367.           if (rsq <= cut_lj_innersq) {
  368.             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  369.            
  370.           } else {
  371.             r3inv = sqrt(r6inv);
  372.  
  373.             /* Steinbach and Brooks (1993) Eqs. 10b and 13 */
  374.             k12 =  lj3[itype][jtype] * cut_ljsq_6 / (cut_ljsq_6-cut_lj_innersq_6);
  375.             k6  = -lj4[itype][jtype] * cut_ljsq_3 / (cut_ljsq_3-cut_lj_innersq_3);
  376.  
  377.             /* Neg. Gradient of Steinbach and Brooks Eq. 10b  */
  378.             forcelj = 0;
  379.             forcelj += 2*6* k12*r6inv*(r6inv-1/cut_ljsq_6);
  380.             forcelj += 2*3* k6 *r3inv*(r3inv-1/cut_ljsq_3);
  381.           }
  382.  
  383.         } else forcelj = 0.0;
  384.  
  385.  
  386.         fpair = (forcecoul + factor_lj*forcelj) * r2inv;
  387.  
  388.         if (rsq > cut_out_on_sq) {
  389.           rsw = (cut_out_off - sqrt(rsq))/cut_out_diff;
  390.           fpair *= rsw*rsw*(-2.0*rsw+3.0);
  391.         }
  392.  
  393.         f[i][0] += delx*fpair;
  394.         f[i][1] += dely*fpair;
  395.         f[i][2] += delz*fpair;
  396.         if (newton_pair || j < nlocal) {
  397.           f[j][0] -= delx*fpair;
  398.           f[j][1] -= dely*fpair;
  399.           f[j][2] -= delz*fpair;
  400.         }
  401.       }
  402.     }
  403.   }
  404. }
  405.  
  406. /* ---------------------------------------------------------------------- */
  407.  
  408. void PairLJCharmmCoulLong::compute_middle()
  409. {
  410.   int i,j,ii,jj,inum,jnum,itype,jtype;
  411.   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
  412.   double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
  413.   double philj,switch1,switch2;
  414.   double rsw;
  415.   int *ilist,*jlist,*numneigh,**firstneigh;
  416.  
  417.   double **x = atom->x;
  418.   double **f = atom->f;
  419.   double *q = atom->q;
  420.   int *type = atom->type;
  421.   int nlocal = atom->nlocal;
  422.   double *special_coul = force->special_coul;
  423.   double *special_lj = force->special_lj;
  424.   int newton_pair = force->newton_pair;
  425.   double qqrd2e = force->qqrd2e;
  426.  
  427.   inum = listmiddle->inum;
  428.   ilist = listmiddle->ilist;
  429.   numneigh = listmiddle->numneigh;
  430.   firstneigh = listmiddle->firstneigh;
  431.  
  432.   double cut_in_off = cut_respa[0];
  433.   double cut_in_on = cut_respa[1];
  434.   double cut_out_on = cut_respa[2];
  435.   double cut_out_off = cut_respa[3];
  436.  
  437.   double cut_in_diff = cut_in_on - cut_in_off;
  438.   double cut_out_diff = cut_out_off - cut_out_on;
  439.   double cut_in_off_sq = cut_in_off*cut_in_off;
  440.   double cut_in_on_sq = cut_in_on*cut_in_on;
  441.   double cut_out_on_sq = cut_out_on*cut_out_on;
  442.   double cut_out_off_sq = cut_out_off*cut_out_off;
  443.  
  444.  
  445.   double cut_ljsq_6;
  446.   double cut_ljsq_3;
  447.   double cut_lj_innersq_6;
  448.   double cut_lj_innersq_3;
  449.   double k12;
  450.   double k6;
  451.   double r3inv;
  452.  
  453.   // only need to compute these once
  454.   cut_ljsq_6 = cut_ljsq*cut_ljsq*cut_ljsq;
  455.   cut_ljsq_3 = sqrt(cut_ljsq_6);
  456.   cut_lj_innersq_6 = cut_lj_innersq*cut_lj_innersq*cut_lj_innersq;
  457.   cut_lj_innersq_3 = sqrt(cut_lj_innersq_6);
  458.  
  459.  
  460.  
  461.   // loop over neighbors of my atoms
  462.  
  463.   for (ii = 0; ii < inum; ii++) {
  464.     i = ilist[ii];
  465.     qtmp = q[i];
  466.     xtmp = x[i][0];
  467.     ytmp = x[i][1];
  468.     ztmp = x[i][2];
  469.     itype = type[i];
  470.  
  471.     jlist = firstneigh[i];
  472.     jnum = numneigh[i];
  473.  
  474.     for (jj = 0; jj < jnum; jj++) {
  475.       j = jlist[jj];
  476.       factor_lj = special_lj[sbmask(j)];
  477.       factor_coul = special_coul[sbmask(j)];
  478.       j &= NEIGHMASK;
  479.  
  480.       delx = xtmp - x[j][0];
  481.       dely = ytmp - x[j][1];
  482.       delz = ztmp - x[j][2];
  483.       rsq = delx*delx + dely*dely + delz*delz;
  484.  
  485.       if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
  486.         r2inv = 1.0/rsq;
  487.         forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
  488.         if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
  489.  
  490.  
  491.         /* FORCE CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  492.         if (rsq <= cut_ljsq) {
  493.             r6inv = r2inv*r2inv*r2inv;
  494.             jtype = type[j];
  495.  
  496.           if (rsq <= cut_lj_innersq) {
  497.             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  498.            
  499.           } else {
  500.             r3inv = sqrt(r6inv);
  501.  
  502.             /* Steinbach and Brooks (1993) Eqs. 10b and 13 */
  503.             k12 =  lj3[itype][jtype] * cut_ljsq_6 / (cut_ljsq_6-cut_lj_innersq_6);
  504.             k6  = -lj4[itype][jtype] * cut_ljsq_3 / (cut_ljsq_3-cut_lj_innersq_3);
  505.  
  506.             /* Neg. Gradient of Steinbach and Brooks Eq. 10b  */
  507.             forcelj = 0;
  508.             forcelj += 2*6* k12*r6inv*(r6inv-1/cut_ljsq_6);
  509.             forcelj += 2*3* k6 *r3inv*(r3inv-1/cut_ljsq_3);
  510.           }
  511.  
  512.         } else forcelj = 0.0;
  513.  
  514.  
  515.         fpair = (forcecoul + factor_lj*forcelj) * r2inv;
  516.         if (rsq < cut_in_on_sq) {
  517.           rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
  518.           fpair *= rsw*rsw*(-2.0*rsw+3.0);
  519.         }
  520.         if (rsq > cut_out_on_sq) {
  521.           rsw = (cut_out_off - sqrt(rsq))/cut_out_diff;
  522.           fpair *= rsw*rsw*(-2.0*rsw+3.0);
  523.         }
  524.  
  525.         f[i][0] += delx*fpair;
  526.         f[i][1] += dely*fpair;
  527.         f[i][2] += delz*fpair;
  528.         if (newton_pair || j < nlocal) {
  529.           f[j][0] -= delx*fpair;
  530.           f[j][1] -= dely*fpair;
  531.           f[j][2] -= delz*fpair;
  532.         }
  533.       }
  534.     }
  535.   }
  536. }
  537.  
  538. /* ---------------------------------------------------------------------- */
  539.  
  540. void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
  541. {
  542.   int i,j,ii,jj,inum,jnum,itype,jtype,itable;
  543.   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
  544.   double fraction,table;
  545.   double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
  546.   double grij,expm2,prefactor,t,erfc;
  547.   double philj,switch1,switch2;
  548.   double rsw;
  549.   int *ilist,*jlist,*numneigh,**firstneigh;
  550.   double rsq;
  551.  
  552.   evdwl = ecoul = 0.0;
  553.   if (eflag || vflag) ev_setup(eflag,vflag);
  554.   else evflag = 0;
  555.  
  556.   double **x = atom->x;
  557.   double **f = atom->f;
  558.   double *q = atom->q;
  559.   int *type = atom->type;
  560.   int nlocal = atom->nlocal;
  561.   double *special_coul = force->special_coul;
  562.   double *special_lj = force->special_lj;
  563.   int newton_pair = force->newton_pair;
  564.   double qqrd2e = force->qqrd2e;
  565.  
  566.   inum = listouter->inum;
  567.   ilist = listouter->ilist;
  568.   numneigh = listouter->numneigh;
  569.   firstneigh = listouter->firstneigh;
  570.  
  571.   double cut_in_off = cut_respa[2];
  572.   double cut_in_on = cut_respa[3];
  573.  
  574.   double cut_in_diff = cut_in_on - cut_in_off;
  575.   double cut_in_off_sq = cut_in_off*cut_in_off;
  576.   double cut_in_on_sq = cut_in_on*cut_in_on;
  577.  
  578.  
  579.   double cut_ljsq_6;
  580.   double cut_ljsq_3;
  581.   double cut_lj_innersq_6;
  582.   double cut_lj_innersq_3;
  583.   double k12;
  584.   double k6;
  585.   double r3inv;
  586.  
  587.   // only need to compute these once
  588.   cut_ljsq_6 = cut_ljsq*cut_ljsq*cut_ljsq;
  589.   cut_ljsq_3 = sqrt(cut_ljsq_6);
  590.   cut_lj_innersq_6 = cut_lj_innersq*cut_lj_innersq*cut_lj_innersq;
  591.   cut_lj_innersq_3 = sqrt(cut_lj_innersq_6);
  592.  
  593.  
  594.  
  595.   // loop over neighbors of my atoms
  596.  
  597.   for (ii = 0; ii < inum; ii++) {
  598.     i = ilist[ii];
  599.     qtmp = q[i];
  600.     xtmp = x[i][0];
  601.     ytmp = x[i][1];
  602.     ztmp = x[i][2];
  603.     itype = type[i];
  604.     jlist = firstneigh[i];
  605.     jnum = numneigh[i];
  606.  
  607.     for (jj = 0; jj < jnum; jj++) {
  608.       j = jlist[jj];
  609.       factor_lj = special_lj[sbmask(j)];
  610.       factor_coul = special_coul[sbmask(j)];
  611.       j &= NEIGHMASK;
  612.  
  613.       delx = xtmp - x[j][0];
  614.       dely = ytmp - x[j][1];
  615.       delz = ztmp - x[j][2];
  616.       rsq = delx*delx + dely*dely + delz*delz;
  617.       jtype = type[j];
  618.  
  619.       if (rsq < cut_bothsq) {
  620.         r2inv = 1.0/rsq;
  621.  
  622.         if (rsq < cut_coulsq) {
  623.           if (!ncoultablebits || rsq <= tabinnersq) {
  624.             r = sqrt(rsq);
  625.             grij = g_ewald * r;
  626.             expm2 = exp(-grij*grij);
  627.             t = 1.0 / (1.0 + EWALD_P*grij);
  628.             erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
  629.             prefactor = qqrd2e * qtmp*q[j]/r;
  630.             forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
  631.             if (rsq > cut_in_off_sq) {
  632.               if (rsq < cut_in_on_sq) {
  633.                 rsw = (r - cut_in_off)/cut_in_diff;
  634.                 forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw);
  635.                 if (factor_coul < 1.0)
  636.                   forcecoul -=
  637.                     (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw);
  638.               } else {
  639.                 forcecoul += prefactor;
  640.                 if (factor_coul < 1.0)
  641.                   forcecoul -= (1.0-factor_coul)*prefactor;
  642.               }
  643.             }
  644.           } else {
  645.             union_int_float_t rsq_lookup;
  646.             rsq_lookup.f = rsq;
  647.             itable = rsq_lookup.i & ncoulmask;
  648.             itable >>= ncoulshiftbits;
  649.             fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
  650.             table = ftable[itable] + fraction*dftable[itable];
  651.             forcecoul = qtmp*q[j] * table;
  652.             if (factor_coul < 1.0) {
  653.               table = ctable[itable] + fraction*dctable[itable];
  654.               prefactor = qtmp*q[j] * table;
  655.               forcecoul -= (1.0-factor_coul)*prefactor;
  656.             }
  657.           }
  658.         } else forcecoul = 0.0;
  659.  
  660.         if (rsq < cut_ljsq && rsq > cut_in_off_sq) {
  661.  
  662.  
  663.           /* FORCE CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  664.           if (rsq <= cut_ljsq) {
  665.             r6inv = r2inv*r2inv*r2inv;
  666.             jtype = type[j];
  667.  
  668.             if (rsq <= cut_lj_innersq) {
  669.               forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  670.            
  671.             } else {
  672.               r3inv = sqrt(r6inv);
  673.  
  674.               /* Steinbach and Brooks (1993) Eqs. 10b and 13 */
  675.               k12 =  lj3[itype][jtype] * cut_ljsq_6 / (cut_ljsq_6-cut_lj_innersq_6);
  676.               k6  = -lj4[itype][jtype] * cut_ljsq_3 / (cut_ljsq_3-cut_lj_innersq_3);
  677.  
  678.               /* Neg. Gradient of Steinbach and Brooks Eq. 10b  */
  679.               forcelj = 0;
  680.               forcelj += 2*6* k12*r6inv*(r6inv-1/cut_ljsq_6);
  681.               forcelj += 2*3* k6 *r3inv*(r3inv-1/cut_ljsq_3);
  682.             }
  683.  
  684.           } else forcelj = 0.0;
  685.  
  686.  
  687.           if (rsq < cut_in_on_sq) {
  688.             rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
  689.             forcelj *= rsw*rsw*(-2.0*rsw+3.0);
  690.           }
  691.         } else forcelj = 0.0;
  692.  
  693.         fpair = (forcecoul + forcelj) * r2inv;
  694.  
  695.         f[i][0] += delx*fpair;
  696.         f[i][1] += dely*fpair;
  697.         f[i][2] += delz*fpair;
  698.         if (newton_pair || j < nlocal) {
  699.           f[j][0] -= delx*fpair;
  700.           f[j][1] -= dely*fpair;
  701.           f[j][2] -= delz*fpair;
  702.         }
  703.  
  704.         if (eflag) {
  705.           if (rsq < cut_coulsq) {
  706.             if (!ncoultablebits || rsq <= tabinnersq) {
  707.               ecoul = prefactor*erfc;
  708.               if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
  709.             } else {
  710.               table = etable[itable] + fraction*detable[itable];
  711.               ecoul = qtmp*q[j] * table;
  712.               if (factor_coul < 1.0) {
  713.                 table = ptable[itable] + fraction*dptable[itable];
  714.                 prefactor = qtmp*q[j] * table;
  715.                 ecoul -= (1.0-factor_coul)*prefactor;
  716.               }
  717.             }
  718.           } else ecoul = 0.0;
  719.  
  720. //          if (rsq < cut_ljsq) {
  721. //            r6inv = r2inv*r2inv*r2inv;
  722. //            evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
  723. //            if (rsq > cut_lj_innersq) {
  724. //              switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  725. //                (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  726. //              evdwl *= switch1;
  727. //            }
  728. //            evdwl *= factor_lj;
  729. //          } else evdwl = 0.0;
  730.  
  731.  
  732.           /* ENERGY CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  733.           if (rsq <= cut_ljsq) {
  734.             if (rsq <= cut_lj_innersq) {
  735.               /* Steinbach and Brooks (1993) Eqs. 10a and 12 */
  736.               philj  = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
  737.               philj += (-lj3[itype][jtype] / (cut_ljsq_6*cut_lj_innersq_6));
  738.               philj += ( lj4[itype][jtype] / (cut_ljsq_3*cut_lj_innersq_3));
  739.             } else {
  740.               /* Steinbach and Brooks (1993) Eqs. 10b */
  741.               philj = 0;
  742.               philj += k12 * (r6inv - 1/cut_ljsq_6) * (r6inv - 1/cut_ljsq_6);
  743.               philj += k6  * (r3inv - 1/cut_ljsq_3) * (r3inv - 1/cut_ljsq_3);
  744.             }
  745.             evdwl = factor_lj*philj;
  746.           } else evdwl = 0.0;
  747.  
  748.  
  749.         }
  750.  
  751.         if (vflag) {
  752.           if (rsq < cut_coulsq) {
  753.             if (!ncoultablebits || rsq <= tabinnersq) {
  754.               forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
  755.               if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
  756.             } else {
  757.               table = vtable[itable] + fraction*dvtable[itable];
  758.               forcecoul = qtmp*q[j] * table;
  759.               if (factor_coul < 1.0) {
  760.                 table = ptable[itable] + fraction*dptable[itable];
  761.                 prefactor = qtmp*q[j] * table;
  762.                 forcecoul -= (1.0-factor_coul)*prefactor;
  763.               }
  764.             }
  765.           } else forcecoul = 0.0;
  766.  
  767.           //~ if (rsq <= cut_in_off_sq) {
  768.             //~ r6inv = r2inv*r2inv*r2inv;
  769.             //~ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  770.             //~ if (rsq > cut_lj_innersq) {
  771.               //~ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  772.                 //~ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  773.               //~ switch2 = 12.0*rsq * (cut_ljsq-rsq) *
  774.                 //~ (rsq-cut_lj_innersq) / denom_lj;
  775.               //~ philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
  776.               //~ forcelj = forcelj*switch1 + philj*switch2;
  777.             //~ }
  778.           //~ } else if (rsq <= cut_in_on_sq) {
  779.             //~ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  780.             //~ if (rsq > cut_lj_innersq) {
  781.               //~ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  782.                 //~ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  783.               //~ switch2 = 12.0*rsq * (cut_ljsq-rsq) *
  784.                 //~ (rsq-cut_lj_innersq) / denom_lj;
  785.               //~ philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
  786.               //~ forcelj = forcelj*switch1 + philj*switch2;
  787.             //~ }
  788.           //~ }
  789.  
  790.           /* FORCE CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  791.           if (rsq <= cut_ljsq) {
  792.             //r6inv = r2inv*r2inv*r2inv;
  793.             //jtype = type[j];
  794.  
  795.             if (rsq <= cut_lj_innersq) {
  796.               forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  797.            
  798.             } else {
  799.               r3inv = sqrt(r6inv);
  800.  
  801.               /* Steinbach and Brooks (1993) Eqs. 10b and 13 */
  802.               k12 =  lj3[itype][jtype] * cut_ljsq_6 / (cut_ljsq_6-cut_lj_innersq_6);
  803.               k6  = -lj4[itype][jtype] * cut_ljsq_3 / (cut_ljsq_3-cut_lj_innersq_3);
  804.  
  805.               /* Neg. Gradient of Steinbach and Brooks Eq. 10b  */
  806.               forcelj = 0;
  807.               forcelj += 2*6* k12*r6inv*(r6inv-1/cut_ljsq_6);
  808.               forcelj += 2*3* k6 *r3inv*(r3inv-1/cut_ljsq_3);
  809.             }
  810.  
  811.           } else forcelj = 0.0;
  812.  
  813.  
  814.  
  815.           fpair = (forcecoul + factor_lj*forcelj) * r2inv;
  816.         }
  817.  
  818.         if (evflag) ev_tally(i,j,nlocal,newton_pair,
  819.                              evdwl,ecoul,fpair,delx,dely,delz);
  820.       }
  821.     }
  822.   }
  823. }
  824.  
  825. /* ----------------------------------------------------------------------
  826.    allocate all arrays
  827. ------------------------------------------------------------------------- */
  828.  
  829. void PairLJCharmmCoulLong::allocate()
  830. {
  831.   allocated = 1;
  832.   int n = atom->ntypes;
  833.  
  834.   memory->create(setflag,n+1,n+1,"pair:setflag");
  835.   for (int i = 1; i <= n; i++)
  836.     for (int j = i; j <= n; j++)
  837.       setflag[i][j] = 0;
  838.  
  839.   memory->create(cutsq,n+1,n+1,"pair:cutsq");
  840.  
  841.   memory->create(epsilon,n+1,n+1,"pair:epsilon");
  842.   memory->create(sigma,n+1,n+1,"pair:sigma");
  843.   memory->create(eps14,n+1,n+1,"pair:eps14");
  844.   memory->create(sigma14,n+1,n+1,"pair:sigma14");
  845.   memory->create(lj1,n+1,n+1,"pair:lj1");
  846.   memory->create(lj2,n+1,n+1,"pair:lj2");
  847.   memory->create(lj3,n+1,n+1,"pair:lj3");
  848.   memory->create(lj4,n+1,n+1,"pair:lj4");
  849.   memory->create(lj14_1,n+1,n+1,"pair:lj14_1");
  850.   memory->create(lj14_2,n+1,n+1,"pair:lj14_2");
  851.   memory->create(lj14_3,n+1,n+1,"pair:lj14_3");
  852.   memory->create(lj14_4,n+1,n+1,"pair:lj14_4");
  853.   memory->create(offset,n+1,n+1,"pair:offset");
  854. }
  855.  
  856. /* ----------------------------------------------------------------------
  857.    global settings
  858.    unlike other pair styles,
  859.      there are no individual pair settings that these override
  860. ------------------------------------------------------------------------- */
  861.  
  862. void PairLJCharmmCoulLong::settings(int narg, char **arg)
  863. {
  864.   if (narg != 2 && narg != 3) error->all(FLERR,"Illegal pair_style command");
  865.  
  866.   cut_lj_inner = force->numeric(FLERR,arg[0]);
  867.   cut_lj = force->numeric(FLERR,arg[1]);
  868.   if (narg == 2) cut_coul = cut_lj;
  869.   else cut_coul = force->numeric(FLERR,arg[2]);
  870. }
  871.  
  872. /* ----------------------------------------------------------------------
  873.    set coeffs for one or more type pairs
  874. ------------------------------------------------------------------------- */
  875.  
  876. void PairLJCharmmCoulLong::coeff(int narg, char **arg)
  877. {
  878.   if (narg != 4 && narg != 6) error->all(FLERR,"Illegal pair_coeff command");
  879.   if (!allocated) allocate();
  880.  
  881.   int ilo,ihi,jlo,jhi;
  882.   force->bounds(arg[0],atom->ntypes,ilo,ihi);
  883.   force->bounds(arg[1],atom->ntypes,jlo,jhi);
  884.  
  885.   double epsilon_one = force->numeric(FLERR,arg[2]);
  886.   double sigma_one = force->numeric(FLERR,arg[3]);
  887.   double eps14_one = epsilon_one;
  888.   double sigma14_one = sigma_one;
  889.   if (narg == 6) {
  890.     eps14_one = force->numeric(FLERR,arg[4]);
  891.     sigma14_one = force->numeric(FLERR,arg[5]);
  892.   }
  893.  
  894.   int count = 0;
  895.   for (int i = ilo; i <= ihi; i++) {
  896.     for (int j = MAX(jlo,i); j <= jhi; j++) {
  897.       epsilon[i][j] = epsilon_one;
  898.       sigma[i][j] = sigma_one;
  899.       eps14[i][j] = eps14_one;
  900.       sigma14[i][j] = sigma14_one;
  901.       setflag[i][j] = 1;
  902.       count++;
  903.     }
  904.   }
  905.  
  906.   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
  907. }
  908.  
  909. /* ----------------------------------------------------------------------
  910.    init specific to this pair style
  911. ------------------------------------------------------------------------- */
  912.  
  913. void PairLJCharmmCoulLong::init_style()
  914. {
  915.   if (!atom->q_flag)
  916.     error->all(FLERR,
  917.                "Pair style lj/charmm/coul/long requires atom attribute q");
  918.  
  919.   // request regular or rRESPA neighbor lists
  920.  
  921.   int irequest;
  922.  
  923.   if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
  924.     int respa = 0;
  925.     if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
  926.     if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
  927.  
  928.     if (respa == 0) irequest = neighbor->request(this,instance_me);
  929.     else if (respa == 1) {
  930.       irequest = neighbor->request(this,instance_me);
  931.       neighbor->requests[irequest]->id = 1;
  932.       neighbor->requests[irequest]->half = 0;
  933.       neighbor->requests[irequest]->respainner = 1;
  934.       irequest = neighbor->request(this,instance_me);
  935.       neighbor->requests[irequest]->id = 3;
  936.       neighbor->requests[irequest]->half = 0;
  937.       neighbor->requests[irequest]->respaouter = 1;
  938.     } else {
  939.       irequest = neighbor->request(this,instance_me);
  940.       neighbor->requests[irequest]->id = 1;
  941.       neighbor->requests[irequest]->half = 0;
  942.       neighbor->requests[irequest]->respainner = 1;
  943.       irequest = neighbor->request(this,instance_me);
  944.       neighbor->requests[irequest]->id = 2;
  945.       neighbor->requests[irequest]->half = 0;
  946.       neighbor->requests[irequest]->respamiddle = 1;
  947.       irequest = neighbor->request(this,instance_me);
  948.       neighbor->requests[irequest]->id = 3;
  949.       neighbor->requests[irequest]->half = 0;
  950.       neighbor->requests[irequest]->respaouter = 1;
  951.     }
  952.  
  953.   } else irequest = neighbor->request(this,instance_me);
  954.  
  955.   // require cut_lj_inner < cut_lj
  956.  
  957.   if (cut_lj_inner >= cut_lj)
  958.     error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
  959.  
  960.   cut_lj_innersq = cut_lj_inner * cut_lj_inner;
  961.   cut_ljsq = cut_lj * cut_lj;
  962.   cut_coulsq = cut_coul * cut_coul;
  963.   cut_bothsq = MAX(cut_ljsq,cut_coulsq);
  964.  
  965.   denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
  966.     (cut_ljsq-cut_lj_innersq);
  967.  
  968.   // set & error check interior rRESPA cutoffs
  969.  
  970.   if (strstr(update->integrate_style,"respa") &&
  971.       ((Respa *) update->integrate)->level_inner >= 0) {
  972.     cut_respa = ((Respa *) update->integrate)->cutoff;
  973.     if (MIN(cut_lj,cut_coul) < cut_respa[3])
  974.       error->all(FLERR,"Pair cutoff < Respa interior cutoff");
  975.     if (cut_lj_inner < cut_respa[1])
  976.       error->all(FLERR,"Pair inner cutoff < Respa interior cutoff");
  977.   } else cut_respa = NULL;
  978.  
  979.   // insure use of KSpace long-range solver, set g_ewald
  980.  
  981.   if (force->kspace == NULL)
  982.     error->all(FLERR,"Pair style requires a KSpace style");
  983.   g_ewald = force->kspace->g_ewald;
  984.  
  985.   // setup force tables
  986.  
  987.   if (ncoultablebits) init_tables(cut_coul,cut_respa);
  988. }
  989.  
  990. /* ----------------------------------------------------------------------
  991.    neighbor callback to inform pair style of neighbor list to use
  992.    regular or rRESPA
  993. ------------------------------------------------------------------------- */
  994.  
  995. void PairLJCharmmCoulLong::init_list(int id, NeighList *ptr)
  996. {
  997.   if (id == 0) list = ptr;
  998.   else if (id == 1) listinner = ptr;
  999.   else if (id == 2) listmiddle = ptr;
  1000.   else if (id == 3) listouter = ptr;
  1001. }
  1002.  
  1003. /* ----------------------------------------------------------------------
  1004.    init for one type pair i,j and corresponding j,i
  1005. ------------------------------------------------------------------------- */
  1006.  
  1007. double PairLJCharmmCoulLong::init_one(int i, int j)
  1008. {
  1009.   if (setflag[i][j] == 0) {
  1010.     epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
  1011.                                sigma[i][i],sigma[j][j]);
  1012.     sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
  1013.     eps14[i][j] = mix_energy(eps14[i][i],eps14[j][j],
  1014.                                sigma14[i][i],sigma14[j][j]);
  1015.     sigma14[i][j] = mix_distance(sigma14[i][i],sigma14[j][j]);
  1016.   }
  1017.  
  1018.   double cut = MAX(cut_lj,cut_coul);
  1019.  
  1020.   lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
  1021.   lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
  1022.   lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
  1023.   lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
  1024.   lj14_1[i][j] = 48.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
  1025.   lj14_2[i][j] = 24.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
  1026.   lj14_3[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
  1027.   lj14_4[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
  1028.  
  1029.   lj1[j][i] = lj1[i][j];
  1030.   lj2[j][i] = lj2[i][j];
  1031.   lj3[j][i] = lj3[i][j];
  1032.   lj4[j][i] = lj4[i][j];
  1033.   lj14_1[j][i] = lj14_1[i][j];
  1034.   lj14_2[j][i] = lj14_2[i][j];
  1035.   lj14_3[j][i] = lj14_3[i][j];
  1036.   lj14_4[j][i] = lj14_4[i][j];
  1037.  
  1038.   return cut;
  1039. }
  1040.  
  1041. /* ----------------------------------------------------------------------
  1042.   proc 0 writes to restart file
  1043. ------------------------------------------------------------------------- */
  1044.  
  1045. void PairLJCharmmCoulLong::write_restart(FILE *fp)
  1046. {
  1047.   write_restart_settings(fp);
  1048.  
  1049.   int i,j;
  1050.   for (i = 1; i <= atom->ntypes; i++)
  1051.     for (j = i; j <= atom->ntypes; j++) {
  1052.       fwrite(&setflag[i][j],sizeof(int),1,fp);
  1053.       if (setflag[i][j]) {
  1054.         fwrite(&epsilon[i][j],sizeof(double),1,fp);
  1055.         fwrite(&sigma[i][j],sizeof(double),1,fp);
  1056.         fwrite(&eps14[i][j],sizeof(double),1,fp);
  1057.         fwrite(&sigma14[i][j],sizeof(double),1,fp);
  1058.       }
  1059.     }
  1060. }
  1061.  
  1062. /* ----------------------------------------------------------------------
  1063.   proc 0 reads from restart file, bcasts
  1064. ------------------------------------------------------------------------- */
  1065.  
  1066. void PairLJCharmmCoulLong::read_restart(FILE *fp)
  1067. {
  1068.   read_restart_settings(fp);
  1069.  
  1070.   allocate();
  1071.  
  1072.   int i,j;
  1073.   int me = comm->me;
  1074.   for (i = 1; i <= atom->ntypes; i++)
  1075.     for (j = i; j <= atom->ntypes; j++) {
  1076.       if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
  1077.       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
  1078.       if (setflag[i][j]) {
  1079.         if (me == 0) {
  1080.           fread(&epsilon[i][j],sizeof(double),1,fp);
  1081.           fread(&sigma[i][j],sizeof(double),1,fp);
  1082.           fread(&eps14[i][j],sizeof(double),1,fp);
  1083.           fread(&sigma14[i][j],sizeof(double),1,fp);
  1084.         }
  1085.         MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
  1086.         MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
  1087.         MPI_Bcast(&eps14[i][j],1,MPI_DOUBLE,0,world);
  1088.         MPI_Bcast(&sigma14[i][j],1,MPI_DOUBLE,0,world);
  1089.       }
  1090.     }
  1091. }
  1092.  
  1093. /* ----------------------------------------------------------------------
  1094.   proc 0 writes to restart file
  1095. ------------------------------------------------------------------------- */
  1096.  
  1097. void PairLJCharmmCoulLong::write_restart_settings(FILE *fp)
  1098. {
  1099.   fwrite(&cut_lj_inner,sizeof(double),1,fp);
  1100.   fwrite(&cut_lj,sizeof(double),1,fp);
  1101.   fwrite(&cut_coul,sizeof(double),1,fp);
  1102.   fwrite(&offset_flag,sizeof(int),1,fp);
  1103.   fwrite(&mix_flag,sizeof(int),1,fp);
  1104.   fwrite(&ncoultablebits,sizeof(int),1,fp);
  1105.   fwrite(&tabinner,sizeof(double),1,fp);
  1106. }
  1107.  
  1108. /* ----------------------------------------------------------------------
  1109.   proc 0 reads from restart file, bcasts
  1110. ------------------------------------------------------------------------- */
  1111.  
  1112. void PairLJCharmmCoulLong::read_restart_settings(FILE *fp)
  1113. {
  1114.   if (comm->me == 0) {
  1115.     fread(&cut_lj_inner,sizeof(double),1,fp);
  1116.     fread(&cut_lj,sizeof(double),1,fp);
  1117.     fread(&cut_coul,sizeof(double),1,fp);
  1118.     fread(&offset_flag,sizeof(int),1,fp);
  1119.     fread(&mix_flag,sizeof(int),1,fp);
  1120.     fread(&ncoultablebits,sizeof(int),1,fp);
  1121.     fread(&tabinner,sizeof(double),1,fp);
  1122.   }
  1123.   MPI_Bcast(&cut_lj_inner,1,MPI_DOUBLE,0,world);
  1124.   MPI_Bcast(&cut_lj,1,MPI_DOUBLE,0,world);
  1125.   MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
  1126.   MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
  1127.   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
  1128.   MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
  1129.   MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
  1130. }
  1131.  
  1132.  
  1133. /* ----------------------------------------------------------------------
  1134.    proc 0 writes to data file
  1135. ------------------------------------------------------------------------- */
  1136.  
  1137. void PairLJCharmmCoulLong::write_data(FILE *fp)
  1138. {
  1139.   for (int i = 1; i <= atom->ntypes; i++)
  1140.     fprintf(fp,"%d %g %g %g %g\n",
  1141.             i,epsilon[i][i],sigma[i][i],eps14[i][i],sigma14[i][i]);
  1142. }
  1143.  
  1144. /* ----------------------------------------------------------------------
  1145.    proc 0 writes all pairs to data file
  1146. ------------------------------------------------------------------------- */
  1147.  
  1148. void PairLJCharmmCoulLong::write_data_all(FILE *fp)
  1149. {
  1150.   for (int i = 1; i <= atom->ntypes; i++)
  1151.     for (int j = i; j <= atom->ntypes; j++)
  1152.       fprintf(fp,"%d %d %g %g %g %g\n",i,j,
  1153.               epsilon[i][j],sigma[i][j],eps14[i][j],sigma14[i][j]);
  1154. }
  1155.  
  1156. /* ---------------------------------------------------------------------- */
  1157.  
  1158. double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
  1159.                                     double rsq,
  1160.                                     double factor_coul, double factor_lj,
  1161.                                     double &fforce)
  1162. {
  1163.   double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
  1164.   double switch1,switch2,fraction,table,forcecoul,forcelj,phicoul,philj;
  1165.   int itable;
  1166.  
  1167.  
  1168.  
  1169.  
  1170.  
  1171.   double cut_ljsq_6;
  1172.   double cut_ljsq_3;
  1173.   double cut_lj_innersq_6;
  1174.   double cut_lj_innersq_3;
  1175.   double k12;
  1176.   double k6;
  1177.   double r3inv;
  1178.  
  1179.   // only need to compute these once
  1180.   cut_ljsq_6 = cut_ljsq*cut_ljsq*cut_ljsq;
  1181.   cut_ljsq_3 = sqrt(cut_ljsq_6);
  1182.   cut_lj_innersq_6 = cut_lj_innersq*cut_lj_innersq*cut_lj_innersq;
  1183.   cut_lj_innersq_3 = sqrt(cut_lj_innersq_6);
  1184.  
  1185.  
  1186.   r2inv = 1.0/rsq;
  1187.   if (rsq < cut_coulsq) {
  1188.     if (!ncoultablebits || rsq <= tabinnersq) {
  1189.       r = sqrt(rsq);
  1190.       grij = g_ewald * r;
  1191.       expm2 = exp(-grij*grij);
  1192.       t = 1.0 / (1.0 + EWALD_P*grij);
  1193.       erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
  1194.       prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
  1195.       forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
  1196.       if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
  1197.     } else {
  1198.       union_int_float_t rsq_lookup;
  1199.       rsq_lookup.f = rsq;
  1200.       itable = rsq_lookup.i & ncoulmask;
  1201.       itable >>= ncoulshiftbits;
  1202.       fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
  1203.       table = ftable[itable] + fraction*dftable[itable];
  1204.       forcecoul = atom->q[i]*atom->q[j] * table;
  1205.       if (factor_coul < 1.0) {
  1206.         table = ctable[itable] + fraction*dctable[itable];
  1207.         prefactor = atom->q[i]*atom->q[j] * table;
  1208.         forcecoul -= (1.0-factor_coul)*prefactor;
  1209.       }
  1210.     }
  1211.   } else forcecoul = 0.0;
  1212.  
  1213. /* EVL
  1214.   if (rsq < cut_ljsq) {
  1215.     r6inv = r2inv*r2inv*r2inv;
  1216.     forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  1217.     if (rsq > cut_lj_innersq) {
  1218.       switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  1219.         (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  1220.       switch2 = 12.0*rsq * (cut_ljsq-rsq) *
  1221.         (rsq-cut_lj_innersq) / denom_lj;
  1222.       philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
  1223.       forcelj = forcelj*switch1 + philj*switch2;
  1224.     }
  1225.   } else forcelj = 0.0;
  1226.   fforce = (forcecoul + factor_lj*forcelj) * r2inv;
  1227. */
  1228.  
  1229.  
  1230.   /* FORCE CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  1231.   if (rsq <= cut_ljsq) {
  1232.     r6inv = r2inv*r2inv*r2inv;
  1233. //    jtype = type[j];
  1234.  
  1235.     if (rsq <= cut_lj_innersq) {
  1236.       forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
  1237.  
  1238.     } else {
  1239.       r3inv = sqrt(r6inv);
  1240.  
  1241.       /* Steinbach and Brooks (1993) Eqs. 10b and 13 */
  1242.       k12 =  lj3[itype][jtype] * cut_ljsq_6 / (cut_ljsq_6-cut_lj_innersq_6);
  1243.       k6  = -lj4[itype][jtype] * cut_ljsq_3 / (cut_ljsq_3-cut_lj_innersq_3);
  1244.  
  1245.       /* Neg. Gradient of Steinbach and Brooks Eq. 10b  */
  1246.       forcelj = 0;
  1247.       forcelj += 2*6* k12*r6inv*(r6inv-1/cut_ljsq_6);
  1248.       forcelj += 2*3* k6 *r3inv*(r3inv-1/cut_ljsq_3);
  1249.     }
  1250.  
  1251.   } else forcelj = 0.0;
  1252.  
  1253.  
  1254.   fforce = (forcecoul + factor_lj*forcelj) * r2inv;
  1255.  
  1256.  
  1257.   double eng = 0.0;
  1258.   if (rsq < cut_coulsq) {
  1259.     if (!ncoultablebits || rsq <= tabinnersq)
  1260.       phicoul = prefactor*erfc;
  1261.     else {
  1262.       table = etable[itable] + fraction*detable[itable];
  1263.       phicoul = atom->q[i]*atom->q[j] * table;
  1264.     }
  1265.     if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
  1266.     eng += phicoul;
  1267.   }
  1268.  
  1269. /* EVL
  1270.   if (rsq < cut_ljsq) {
  1271.     philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
  1272.     if (rsq > cut_lj_innersq) {
  1273.       switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
  1274.         (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
  1275.       philj *= switch1;
  1276.     }
  1277.     eng += factor_lj*philj;
  1278.   }
  1279. */
  1280.  
  1281.  
  1282.   /* ENERGY CALC USING FORCE SWITCHING (FSW2) FOR CHARMM */
  1283.   if (rsq <= cut_ljsq) {
  1284.     if (rsq <= cut_lj_innersq) {
  1285.       /* Steinbach and Brooks (1993) Eqs. 10a and 12 */
  1286.       philj  = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
  1287.       philj += (-lj3[itype][jtype] / (cut_ljsq_6*cut_lj_innersq_6));
  1288.       philj += ( lj4[itype][jtype] / (cut_ljsq_3*cut_lj_innersq_3));
  1289.     } else {
  1290.       /* Steinbach and Brooks (1993) Eqs. 10b */
  1291.       philj = 0;
  1292.       philj += k12 * (r6inv - 1/cut_ljsq_6) * (r6inv - 1/cut_ljsq_6);
  1293.       philj += k6  * (r3inv - 1/cut_ljsq_3) * (r3inv - 1/cut_ljsq_3);
  1294.     }
  1295.     eng += factor_lj*philj;
  1296.   }
  1297.  
  1298.  
  1299.   return eng;
  1300. }
  1301.  
  1302. /* ---------------------------------------------------------------------- */
  1303.  
  1304. void *PairLJCharmmCoulLong::extract(const char *str, int &dim)
  1305. {
  1306.   dim = 2;
  1307.   if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1;
  1308.   if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
  1309.   if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
  1310.   if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
  1311.  
  1312.   dim = 0;
  1313.   if (strcmp(str,"implicit") == 0) return (void *) &implicit;
  1314.   if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
  1315.  
  1316.   return NULL;
  1317. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!