git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4433 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/force.cpp

@@ -48,7 +48,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;
+  special_angle = special_dihedral = 0;
   special_extra = 0;
 
   dielectric = 1.0;
@@ -415,6 +415,12 @@ void Force::set_special(int narg, char **arg)
       special_coul[2] = atof(arg[iarg+2]);
       special_coul[3] = atof(arg[iarg+3]);
       iarg += 4;
+    } else if (strcmp(arg[iarg],"angle") == 0) {
+      if (iarg+2 > narg) error->all("Illegal special_bonds command");
+      if (strcmp(arg[iarg+1],"no") == 0) special_angle = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) special_angle = 1;
+      else error->all("Illegal special_bonds command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"dihedral") == 0) {
       if (iarg+2 > narg) error->all("Illegal special_bonds command");
       if (strcmp(arg[iarg+1],"no") == 0) special_dihedral = 0;

src/force.h

@@ -53,6 +53,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_angle;         // 0 if defined angles are ignored
+                             // 1 if only weight 1,3 atoms if in an angle
   int special_dihedral;      // 0 if defined dihedrals are ignored
                              // 1 if only weight 1,4 atoms if in a dihedral
   int special_extra;         // extra space for added bonds

src/input.cpp

@@ -1140,6 +1140,7 @@ void Input::special_bonds()
   double lj3 = force->special_lj[3];
   double coul2 = force->special_coul[2];
   double coul3 = force->special_coul[3];
+  int angle = force->special_angle;
   int dihedral = force->special_dihedral;
   int extra = force->special_extra;
 
@@ -1150,7 +1151,9 @@ void Input::special_bonds()
   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] ||
-	dihedral != force->special_dihedral || extra != force->special_extra) {
+	angle != force->special_angle || 
+	dihedral != force->special_dihedral || 
+	extra != force->special_extra) {
       Special special(lmp);
       special.build();
     }

src/special.cpp

@@ -387,6 +387,7 @@ void Special::build()
 
   if (force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) {
     combine();
+    if (force->special_angle) angle_trim();
     return;
   }
 
@@ -546,6 +547,7 @@ void Special::build()
   delete [] bufcopy;
 
   combine();
+  if (force->special_angle) angle_trim();
   if (force->special_dihedral) dihedral_trim();
 }
 
@@ -694,6 +696,171 @@ void Special::combine()
   atom->map_set();
 }
 
+/* ----------------------------------------------------------------------
+   trim list of 1-3 neighbors by checking defined angles
+   delete a 1-3 neigh if they are not end atoms of a defined angle
+------------------------------------------------------------------------- */
+
+void Special::angle_trim()
+{
+  int i,j,m,n,iglobal,jglobal,ilocal,jlocal;
+  MPI_Request request;
+  MPI_Status status;
+  
+  int *num_angle = atom->num_angle;
+  int **angle_atom1 = atom->angle_atom1;
+  int **angle_atom3 = atom->angle_atom3;
+  int **nspecial = atom->nspecial;
+  int **special = atom->special;
+  int nlocal = atom->nlocal;
+
+  // stats on old 1-3 neighbor counts
+
+  double onethreecount = 0.0;
+  for (i = 0; i < nlocal; i++)
+    onethreecount += nspecial[i][1] - nspecial[i][0];
+  
+  double allcount;
+  MPI_Allreduce(&onethreecount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
+  
+  if (me == 0) {
+    if (screen)
+      fprintf(screen,
+	      "  %g = # of 1-3 neighbors before angle trim\n",allcount);
+    if (logfile)
+      fprintf(logfile,
+	      "  %g = # of 1-3 neighbors before angle trim\n",allcount);
+  }
+  
+  // if angles are defined, flag each 1-3 neigh if it appears in an angle
+
+  if (num_angle && atom->nangles) {
+
+    // dflag = flag for 1-3 neighs of all owned atoms
+
+    int maxcount = 0;
+    for (i = 0; i < nlocal; i++)
+      maxcount = MAX(maxcount,nspecial[i][1]-nspecial[i][0]);
+    int **dflag =
+      memory->create_2d_int_array(nlocal,maxcount,"special::dflag");
+    
+    for (i = 0; i < nlocal; i++) {
+      n = nspecial[i][1] - nspecial[i][0];
+      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,3 atoms in each angle stored by atom
+    
+    int nbuf = 0;
+    for (i = 0; i < nlocal; i++) nbuf += 2*num_angle[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,3 atoms in each angle
+
+    int size = 0;
+    for (i = 0; i < nlocal; i++)
+      for (j = 0; j < num_angle[i]; j++) {
+	buf[size++] = angle_atom1[i][j];
+	buf[size++] = angle_atom3[i][j];
+      }
+
+    // cycle buffer around ring of procs back to self
+    // when receive buffer, scan list of 1,3 atoms looking for atoms I own
+    // when find one, scan its 1-3 neigh list and mark I,J as in an angle
+
+    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][0]; m < nspecial[ilocal][1]; m++)
+	    if (jglobal == special[ilocal][m]) {
+	      dflag[ilocal][m-nspecial[ilocal][0]] = 1;
+	      break;
+	    }
+	if (jlocal >= 0 && jlocal < nlocal)
+	  for (m = nspecial[jlocal][0]; m < nspecial[jlocal][1]; m++)
+	    if (iglobal == special[jlocal][m]) {
+	      dflag[jlocal][m-nspecial[jlocal][0]] = 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-3 neighbors if they are not flagged in dflag
+    // preserve 1-4 neighbors
+    
+    int offset;
+    for (i = 0; i < nlocal; i++) {
+      offset = m = nspecial[i][0];
+      for (j = nspecial[i][0]; j < nspecial[i][1]; j++)
+	if (dflag[i][j-offset]) special[i][m++] = special[i][j];
+      offset = m;
+      for (j = nspecial[i][1]; j < nspecial[i][2]; j++)
+	special[i][m++] = special[i][j];
+      nspecial[i][1] = offset;
+      nspecial[i][2] = m;
+    }
+    
+    // clean up
+
+    memory->destroy_2d_int_array(dflag);
+    delete [] buf;
+    delete [] bufcopy;
+
+  // if no angles are defined, delete all 1-3 neighs, preserving 1-4 neighs
+
+  } else {
+    for (i = 0; i < nlocal; i++) {
+      m = nspecial[i][0];
+      for (j = nspecial[i][1]; j < nspecial[i][2]; j++)
+	special[i][m++] = special[i][j];
+      nspecial[i][1] = nspecial[i][0];
+      nspecial[i][2] = m;
+    }
+  }
+
+  // stats on new 1-3 neighbor counts
+
+  onethreecount = 0.0;
+  for (i = 0; i < nlocal; i++)
+    onethreecount += nspecial[i][1] - nspecial[i][0];
+  
+  MPI_Allreduce(&onethreecount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
+
+  if (me == 0) {
+    if (screen)
+      fprintf(screen,
+	      "  %g = # of 1-3 neighbors after angle trim\n",allcount);
+    if (logfile)
+      fprintf(logfile,
+	      "  %g = # of 1-3 neighbors after angle trim\n",allcount);
+  }  
+}
+
 /* ----------------------------------------------------------------------
    trim list of 1-4 neighbors by checking defined dihedrals
    delete a 1-4 neigh if they are not end atoms of a defined dihedral

src/special.h

@@ -30,6 +30,7 @@ class Special : protected Pointers {
   int dihedral_flag;
 
   void combine();
+  void angle_trim();
   void dihedral_trim();
 };