Commit 4baa665a authored by lucienPan's avatar lucienPan
Browse files

Patch of class2 dihedral 

Fix the NAN problem when any two bonds are nearly parallel
parent 90729ebe

src/CLASS2/dihedral_class2.cpp

+44 −40
Original line number Diff line number Diff line
@@ -27,7 +27,6 @@ 
#include "math_const.h"

#include "memory.h"

#include "error.h"

#include "utils.h"



using namespace LAMMPS_NS;

using namespace MathConst;
@@ -182,6 +181,11 @@ void DihedralClass2::compute(int eflag, int vflag) 
    costh13 = c0;

    costh23 = (vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z) * r12c2;



    costh12 = MAX(MIN(costh12, 1.0), -1.0);

    costh13 = MAX(MIN(costh13, 1.0), -1.0);

    costh23 = MAX(MIN(costh23, 1.0), -1.0);

    c0 = costh13;



    // cos and sin of 2 angles and final c



    sin2 = MAX(1.0 - costh12*costh12,0.0);
@@ -836,45 +840,45 @@ void DihedralClass2::read_restart(FILE *fp) 
  allocate();



  if (comm->me == 0) {

    utils::sfread(FLERR,&k1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&k2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&k3[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&phi1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&phi2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&phi3[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&mbt_f1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&mbt_f2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&mbt_f3[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&mbt_r0[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&ebt_f1_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&ebt_f2_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&ebt_f3_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&ebt_r0_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&ebt_f1_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&ebt_f2_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&ebt_f3_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&ebt_r0_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&at_f1_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&at_f2_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&at_f3_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&at_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&at_f1_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&at_f2_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&at_f3_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&at_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&aat_k[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&aat_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&aat_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);



    utils::sfread(FLERR,&bb13t_k[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&bb13t_r10[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    utils::sfread(FLERR,&bb13t_r30[1],sizeof(double),atom->ndihedraltypes,fp,NULL,error);

    fread(&k1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&k2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&k3[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&phi1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&phi2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&phi3[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&mbt_f1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&mbt_f2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&mbt_f3[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&mbt_r0[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&ebt_f1_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&ebt_f2_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&ebt_f3_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&ebt_r0_1[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&ebt_f1_2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&ebt_f2_2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&ebt_f3_2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&ebt_r0_2[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&at_f1_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&at_f2_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&at_f3_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&at_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&at_f1_2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&at_f2_2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&at_f3_2[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&at_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&aat_k[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&aat_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&aat_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp);



    fread(&bb13t_k[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&bb13t_r10[1],sizeof(double),atom->ndihedraltypes,fp);

    fread(&bb13t_r30[1],sizeof(double),atom->ndihedraltypes,fp);

  }



  MPI_Bcast(&k1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);

CRA Git | Maintained and supported by SUSTech CRA and CCSE