Support CHARMM psf that contains colinear type of lonepair.

wrappers/python/simtk/openmm/app/charmmpsffile.py

@@ -366,25 +366,30 @@ class CharmmPsfFile(object):
         holder = psfsections['NUMLP NUMLPH'][1]
         lonepair_list = TrackedList()
         if numlp != 0 or numlph != 0:
+            lp_hostnum_list=[]
             lp_distance_list=[]
             lp_angle_list=[]
             lp_dihe_list=[]
             for i in range(numlp):
                 lpline = holder[i].split()
-                if len(lpline)!=6 or lpline[0] != '3' or lpline[2] != 'F' or int(lpline[1]) != 4*i+1 :
+                if len(lpline)!=6 or lpline[2] != 'F' :
                     raise CharmmPSFError('Lonepair format error')
                 else:
+                    lp_hostnum_list.append(int(lpline[0]))
                     lp_distance_list.append(float(lpline[3]))
                     lp_angle_list.append(float(lpline[4]))
                     lp_dihe_list.append(float(lpline[5]))
+            lp_atom_counter=0
             for i in range(numlp):
-                lpline = holder[int(i/2)+numlp].split()
-                icolumn = (i%2) * 4
-                id1=int(lpline[icolumn])  -1
-                id2=int(lpline[icolumn+1])-1
-                id3=int(lpline[icolumn+2])-1
-                id4=int(lpline[icolumn+3])-1
-                lonepair_list.append([id1, id2, id3, id4, lp_distance_list[i], lp_angle_list[i], lp_dihe_list[i]])
+                idall=[]
+                for j in range(lp_hostnum_list[i]+1):
+                    iline = int((lp_atom_counter+j)/8)+numlp
+                    icolumn = (lp_atom_counter+j)%8
+                    idall.append(int(holder[iline].split()[icolumn])-1)
+                if len(idall)==3:
+                    idall.append(-1) # use id4=-1 to mark colinear
+                lonepair_list.append([idall[0], idall[1], idall[2], idall[3], lp_distance_list[i], lp_angle_list[i], lp_dihe_list[i]])
+                lp_atom_counter += lp_hostnum_list[i]+1
         # In Drude psf, here comes anisotropic section
         if IsDrudePSF:
             numaniso = psfsections['NUMANISO'][0]
@@ -912,17 +917,27 @@ class CharmmPsfFile(object):
             atom1=lpsite[1]
             atom2=lpsite[2]
             atom3=lpsite[3]
-            if lpsite[4] > 0 : #relative
-                r = lpsite[4] /10.0 # in nanometer
-                xweights = [-1.0, 0.0, 1.0]
-            elif lpsite[4] < 0: # bisector
-                r = lpsite[4] / (-10.0)
-                xweights = [-1.0, 0.5, 0.5]
-            theta = lpsite[5] * pi / 180.0
-            phi = (180.0 - lpsite[6]) * pi / 180.0
-            p = [r*cos(theta), r*sin(theta)*cos(phi), r*sin(theta)*sin(phi)]
-            p = [x if abs(x) > 1e-10 else 0 for x in p] # Avoid tiny numbers caused by roundoff error
-            system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom3, atom2, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(xweights[0],xweights[1],xweights[2]), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(p[0],p[1],p[2])))
+            if atom3 > 0: 
+                if lpsite[4] > 0 : # relative lonepair type
+                    r = lpsite[4] /10.0 # in nanometer
+                    xweights = [-1.0, 0.0, 1.0]
+                elif lpsite[4] < 0: # bisector lonepair type
+                    r = lpsite[4] / (-10.0)
+                    xweights = [-1.0, 0.5, 0.5]
+                theta = lpsite[5] * pi / 180.0
+                phi = (180.0 - lpsite[6]) * pi / 180.0
+                p = [r*cos(theta), r*sin(theta)*cos(phi), r*sin(theta)*sin(phi)]
+                p = [x if abs(x) > 1e-10 else 0 for x in p] # Avoid tiny numbers caused by roundoff error
+                system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom3, atom2, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(xweights[0],xweights[1],xweights[2]), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(p[0],p[1],p[2])))
+            else: # colinear lonepair type
+                # find a real atom to be the third one for LocalCoordinatesSite
+                for bond in self.bond_list:
+                    if (bond.atom1.idx == atom2 and bond.atom2.idx != atom1):
+                        a3=bond.atom2.idx
+                    elif (bond.atom2.idx == atom2 and bond.atom1.idx != atom1):
+                        a3=bond.atom1.idx
+                r = lpsite[4] / 10.0 # in nanometer
+                system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom2, a3, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(1.0,-1.0, 0.0), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(r,0.0,0.0)))   
         # Add Bond forces
         if verbose: print('Adding bonds...')
         force = mm.HarmonicBondForce()