From 71fbecea76ac03565ffbfd36352f3869d8460123 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 4 Dec 2008 19:47:29 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2286 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/force.cpp | 49 +++++++++++---- src/force.h | 2 + src/input.cpp | 8 ++- src/special.cpp | 163 ++++++++++++++++++++++++++++++++++++++++++++++++ src/special.h | 2 + 5 files changed, 208 insertions(+), 16 deletions(-) diff --git a/src/force.cpp b/src/force.cpp index f394db452c..400840b143 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -55,6 +55,7 @@ Force::Force(LAMMPS *lmp) : Pointers(lmp) special_lj[1] = special_lj[2] = special_lj[3] = 0.0; special_coul[1] = special_coul[2] = special_coul[3] = 0.0; + special_dihedral = 0; dielectric = 1.0; @@ -352,31 +353,53 @@ void Force::set_special(int narg, char **arg) { if (narg == 0) error->all("Illegal special_bonds command"); - if (narg == 1 && strcmp(arg[0],"charmm") == 0) { + if (strcmp(arg[0],"charmm") == 0) { + if (narg != 1) error->all("Illegal special_bonds command"); special_lj[1] = 0.0; special_lj[2] = 0.0; special_lj[3] = 0.0; special_coul[1] = 0.0; special_coul[2] = 0.0; special_coul[3] = 0.0; - } else if (narg == 1 && strcmp(arg[0],"amber") == 0) { + special_dihedral = 0; + return; + } + + if (strcmp(arg[0],"amber") == 0) { + if (narg != 1) error->all("Illegal special_bonds command"); special_lj[1] = 0.0; special_lj[2] = 0.0; special_lj[3] = 0.5; special_coul[1] = 0.0; special_coul[2] = 0.0; special_coul[3] = 5.0/6.0; - } else if (narg == 3) { - special_lj[1] = special_coul[1] = atof(arg[0]); - special_lj[2] = special_coul[2] = atof(arg[1]); - special_lj[3] = special_coul[3] = atof(arg[2]); - } else if (narg == 6) { - special_lj[1] = atof(arg[0]); - special_lj[2] = atof(arg[1]); - special_lj[3] = atof(arg[2]); - special_coul[1] = atof(arg[3]); - special_coul[2] = atof(arg[4]); - special_coul[3] = atof(arg[5]); + special_dihedral = 0; + return; + } + + int iarg; + if (strcmp(arg[0],"dihedral") == 0) { + special_dihedral = 1; + iarg = 1; + } else if (strcmp(arg[0],"explicit") == 0) { + special_dihedral = 0; + iarg = 1; + } else { + special_dihedral = 0; + iarg = 0; + } + + if (narg-iarg == 3) { + special_lj[1] = special_coul[1] = atof(arg[iarg+0]); + special_lj[2] = special_coul[2] = atof(arg[iarg+1]); + special_lj[3] = special_coul[3] = atof(arg[iarg+2]); + } else if (narg-iarg == 6) { + special_lj[1] = atof(arg[iarg+0]); + special_lj[2] = atof(arg[iarg+1]); + special_lj[3] = atof(arg[iarg+2]); + special_coul[1] = atof(arg[iarg+3]); + special_coul[2] = atof(arg[iarg+4]); + special_coul[3] = atof(arg[iarg+5]); } else error->all("Illegal special_bonds command"); } diff --git a/src/force.h b/src/force.h index 6bbb566488..5cc62080d0 100644 --- a/src/force.h +++ b/src/force.h @@ -52,6 +52,8 @@ class Force : protected Pointers { // index [0] is not used in these arrays double special_lj[4]; // 1-2, 1-3, 1-4 prefactors for LJ double special_coul[4]; // 1-2, 1-3, 1-4 prefactors for Coulombics + int special_dihedral; // 0 if defined dihedrals are ignored + // 1 if only weight 1,4 atoms if in a dihedral Force(class LAMMPS *); ~Force(); diff --git a/src/input.cpp b/src/input.cpp index 54acb3aba4..4872e86d68 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -1080,20 +1080,22 @@ void Input::shape() void Input::special_bonds() { - // store 1-3,1-4 values before change + // store 1-3,1-4 and dihedral flag values before change double lj2 = force->special_lj[2]; double lj3 = force->special_lj[3]; double coul2 = force->special_coul[2]; double coul3 = force->special_coul[3]; + int dihedral = force->special_dihedral; force->set_special(narg,arg); - // if simulation box defined and 1-3,1-4 values changed, redo special list + // if simulation box defined and saved values changed, redo special list if (domain->box_exist && atom->molecular) { if (lj2 != force->special_lj[2] || lj3 != force->special_lj[3] || - coul2 != force->special_coul[2] || coul3 != force->special_coul[3]) { + coul2 != force->special_coul[2] || coul3 != force->special_coul[3] || + dihedral != force->special_dihedral) { Special special(lmp); special.build(); } diff --git a/src/special.cpp b/src/special.cpp index c1ddbfbe88..7b66f60ad0 100644 --- a/src/special.cpp +++ b/src/special.cpp @@ -546,6 +546,7 @@ void Special::build() delete [] bufcopy; combine(); + if (force->special_dihedral) dihedral_trim(); } /* ---------------------------------------------------------------------- @@ -690,3 +691,165 @@ void Special::combine() atom->nghost = 0; atom->map_set(); } + +/* ---------------------------------------------------------------------- + trim list of 1-4 neighbors by checking defined dihedrals + delete a 1-4 neigh if those atoms are not end points of a defined dihedral +------------------------------------------------------------------------- */ + +void Special::dihedral_trim() +{ + int i,j,m,n,iglobal,jglobal,ilocal,jlocal; + MPI_Request request; + MPI_Status status; + + int *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + int **dihedral_atom1 = atom->dihedral_atom1; + int **dihedral_atom4 = atom->dihedral_atom4; + int **nspecial = atom->nspecial; + int **special = atom->special; + int nlocal = atom->nlocal; + + // stats on old 1-4 neighbor counts + + double onefourcount = 0.0; + for (i = 0; i < nlocal; i++) + onefourcount += nspecial[i][2] - nspecial[i][1]; + + double allcount; + MPI_Allreduce(&onefourcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world); + + if (me == 0) { + if (screen) + fprintf(screen, + " %g = # of 1-4 neighbors before dihedral trim\n",allcount); + if (logfile) + fprintf(logfile, + " %g = # of 1-4 neighbors before dihedral trim\n",allcount); + } + + // if dihedrals are defined, flag each 1-4 neigh if it appears in a dihedral + + if (num_dihedral && atom->ndihedrals) { + + // dflag = flag for 1-4 neighs of all owned atoms + + int maxcount = 0; + for (i = 0; i < nlocal; i++) + maxcount = MAX(maxcount,nspecial[i][2]-nspecial[i][1]); + int **dflag = + memory->create_2d_int_array(nlocal,maxcount,"special::dflag"); + + for (i = 0; i < nlocal; i++) { + n = nspecial[i][2] - nspecial[i][1]; + for (j = 0; j < n; j++) dflag[i][j] = 0; + } + + // nbufmax = largest buffer needed to hold info from any proc + // info for each atom = list of 1,4 atoms in each dihedral stored by atom + + int nbuf = 0; + for (i = 0; i < nlocal; i++) nbuf += 2*num_dihedral[i]; + + int nbufmax; + MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world); + + int *buf = new int[nbufmax]; + int *bufcopy = new int[nbufmax]; + + // fill buffer with list of 1,4 atoms in each dihedral + + int size = 0; + for (i = 0; i < nlocal; i++) + for (j = 0; j < num_dihedral[i]; j++) { + buf[size++] = dihedral_atom1[i][j]; + buf[size++] = dihedral_atom4[i][j]; + } + + // cycle buffer around ring of procs back to self + // when receive buffer, scan list of 1,4 atoms looking for atoms I own + // when find one, scan its 1-4 neigh list and mark I,J as in a dihedral + + int next = me + 1; + int prev = me -1; + if (next == nprocs) next = 0; + if (prev < 0) prev = nprocs - 1; + + int messtag = 7; + for (int loop = 0; loop < nprocs; loop++) { + i = 0; + while (i < size) { + iglobal = buf[i]; + jglobal = buf[i+1]; + ilocal = atom->map(iglobal); + jlocal = atom->map(jglobal); + if (ilocal >= 0 && ilocal < nlocal) + for (m = nspecial[ilocal][1]; m < nspecial[ilocal][2]; m++) + if (jglobal == special[ilocal][m]) { + dflag[ilocal][m-nspecial[ilocal][1]] = 1; + break; + } + if (jlocal >= 0 && jlocal < nlocal) + for (m = nspecial[jlocal][1]; m < nspecial[jlocal][2]; m++) + if (iglobal == special[jlocal][m]) { + dflag[jlocal][m-nspecial[jlocal][1]] = 1; + break; + } + i += 2; + } + if (me != next) { + MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request); + MPI_Send(buf,size,MPI_INT,next,messtag,world); + MPI_Wait(&request,&status); + MPI_Get_count(&status,MPI_INT,&size); + for (j = 0; j < size; j++) buf[j] = bufcopy[j]; + } + } + + // delete 1-4 neighbors if they are not flagged in dflag + + int offset; + for (i = 0; i < nlocal; i++) { + offset = m = nspecial[i][1]; + for (j = nspecial[i][1]; j < nspecial[i][2]; j++) + if (dflag[i][j-offset]) special[i][m++] = special[i][j]; + nspecial[i][2] = m; + } + + // clean up + + memory->destroy_2d_int_array(dflag); + delete [] buf; + delete [] bufcopy; + + // if no dihedrals are defined, delete all 1-4 neighs + + } else for (i = 0; i < nlocal; i++) nspecial[i][2] = nspecial[i][1]; + + // reset atom->maxspecial + + int maxspecial = 0; + for (i = 0; i < nlocal; i++) + maxspecial = MAX(maxspecial,nspecial[i][2]); + + MPI_Allreduce(&maxspecial,&atom->maxspecial,1,MPI_INT,MPI_MAX,world); + atom->maxspecial = MAX(atom->maxspecial,1); + + // stats on new 1-4 neighbor counts + + onefourcount = 0.0; + for (i = 0; i < nlocal; i++) + onefourcount += nspecial[i][2] - nspecial[i][1]; + + MPI_Allreduce(&onefourcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world); + + if (me == 0) { + if (screen) + fprintf(screen, + " %g = # of 1-4 neighbors after dihedral trim\n",allcount); + if (logfile) + fprintf(logfile, + " %g = # of 1-4 neighbors after dihedral trim\n",allcount); + } +} diff --git a/src/special.h b/src/special.h index 96a03f8aad..fbae488eb7 100644 --- a/src/special.h +++ b/src/special.h @@ -27,8 +27,10 @@ class Special : protected Pointers { private: int me,nprocs; int **onetwo,**onethree,**onefour; + int dihedral_flag; void combine(); + void dihedral_trim(); }; }