1 subroutine pre_collaps(xs,id,ndx,vec,msa,nsa,touch,mode,nmoves, 2 + ncolgr,icolgr,lfnout,lfnpdb,lrgpdb,sysnam,filpdb,iopt,box, 3 + num,amass,mat, 4 + csa,isat,isgm,imol,ifra,vs, 5 + cwa,iwat,xw,vw,mwm,mwa,nwm,nwa, 6 + xwc,vwc,mwmc,nwmc,slvnam,iropt,irrand,nxrep,nyrep,nzrep,drep, 7 + msb,nsb,idsb,rdist,nskip,iskip,lang,lfnmrg,nmerge,xmerge,filmrg, 8 + irenum,invert,ihop,ips) 9c 10 implicit none 11c 12#include "mafdecls.fh" 13 logical pre_wrtpdb 14 external pre_wrtpdb 15c 16 integer msa,nsa,mode,nmoves 17 integer id(msa),ndx(msa,3) 18 real*8 xs(3,msa),vec(6,msa) 19 real*8 touch 20 integer ncolgr 21 integer icolgr(2,ncolgr) 22c 23 integer lfnout,lfnpdb,iopt,mwm,nwm,mwa,nwa,mat,mwmc,nwmc 24 integer lfnmrg,lrgpdb 25 character*255 filpdb,filmrg(100) 26 integer num(mat),isat(msa),isgm(msa),imol(msa),ifra(msa),iwat(mwa) 27 integer ihop(msa),ips(msa) 28 character*16 cwa(mwa),csa(msa) 29 character*3 slvnam 30 real*8 amass(mat),vs(3,msa),xw(3,mwa,mwm),vw(3,mwa,mwm) 31 real*8 xwc(3,mwa,mwmc),vwc(3,mwa,mwmc),box(3) 32 integer iropt,irrand,nxrep,nyrep,nzrep,nskip 33 real*8 rdist,drep(3),xmerge(3,100) 34 integer msb,nsb,nmols,irenum,invert,nmerge 35 integer idsb(2,msb),iskip(3,nskip),lang(*) 36 character*80 sysnam 37c 38 character*255 filnam 39c 40 logical first,found 41 integer i,j,k,m 42 integer nmol,large,nummov 43 real*8 rt,rk,dist 44 real*8 a,b,c,d,e,r,r1,r2,root,dt,dtran 45 real*8 tr(3),tran(3),xmax(3),xmin(3) 46 real*8 a2,b2,c2,d2,p(3),factor 47c 48 write(*,'(/,a,i4,f12.6,/)') ' Collapsing using mode ',mode,touch 49c 50c find total number of molecules 51c 52c on entry, id(i) contains the molecule number of atom i 53c 54 do 13 i=1,nsa 55 ndx(i,1)=0 56 ndx(i,2)=0 57 ndx(i,3)=0 58 13 continue 59c 60c store in ndx(i,1) the number of atoms in molecule i 61c 62 nmol=0 63 do 1 i=1,nsa 64 nmol=max(nmol,id(i)) 65 ndx(id(i),1)=ndx(id(i),1)+1 66 1 continue 67c 68 do 14 i=1,nsa 69 ndx(i,2)=id(i) 70 14 continue 71c 72c group single atom molecules with nearest multi-atom molecules 73c 74 do 15 i=1,nsa 75 if(ndx(id(i),1).eq.1) then 76 k=0 77 rk=1.0d10 78 do 16 j=1,nsa 79 if(ndx(id(j),1).gt.1) then 80 dist=(xs(1,i)-xs(1,j))**2+(xs(2,i)-xs(2,j))**2+ 81 + (xs(3,i)-xs(3,j))**2 82 if(dist.lt.rk) then 83 k=j 84 rk=dist 85 endif 86 endif 87 16 continue 88 if(k.gt.0) then 89 ndx(id(i),1)=0 90 ndx(id(k),1)=ndx(id(k),1)+1 91 ndx(i,2)=id(k) 92 write(*,'(a,i8,a,i4,a,f12.6)') 93 + ' Single ion ',i,' grouped with molecule ', 94 + id(k),', distance is ',sqrt(rk) 95 endif 96 endif 97 15 continue 98c 99 if(ncolgr.gt.0) then 100 do 26 j=1,ncolgr 101 do 27 i=1,nsa 102 if(ndx(i,2).eq.icolgr(2,j)) ndx(i,2)=icolgr(1,j) 103 27 continue 104 write(*,'(a,i8,a,i4)') 105 + ' Molecule ',icolgr(2,j), 106 + ' grouped with molecule ',icolgr(1,j) 107 26 continue 108 endif 109c 110 nummov=0 111c 112 2 continue 113c 114c find molecule with largest possible tranlation 115c 116 large=0 117 dtran=0.0d0 118 root=0.0d0 119c 120c determine centers of geometry 121c 122 do 31 m=1,nmol 123 vec(1,m)=0.0d0 124 vec(2,m)=0.0d0 125 vec(3,m)=0.0d0 126 vec(4,m)=0.0d0 127 vec(5,m)=0.0d0 128 vec(6,m)=0.0d0 129 31 continue 130 do 32 i=1,nsa 131 m=ndx(i,2) 132 vec(1,m)=xs(1,i) 133 vec(2,m)=xs(2,i) 134 vec(3,m)=xs(3,i) 135 vec(4,m)=xs(1,i) 136 vec(5,m)=xs(2,i) 137 vec(6,m)=xs(3,i) 138 32 continue 139 do 33 i=1,nsa 140 m=ndx(i,2) 141 vec(1,m)=min(vec(1,m),xs(1,i)) 142 vec(2,m)=min(vec(2,m),xs(2,i)) 143 vec(3,m)=min(vec(3,m),xs(3,i)) 144 vec(4,m)=max(vec(4,m),xs(1,i)) 145 vec(5,m)=max(vec(5,m),xs(2,i)) 146 vec(6,m)=max(vec(6,m),xs(3,i)) 147 33 continue 148 do 34 m=1,nmol 149 if(ndx(m,1).gt.0) then 150 vec(1,m)=0.5d0*(vec(1,m)+vec(4,m)) 151 vec(2,m)=0.5d0*(vec(2,m)+vec(5,m)) 152 vec(3,m)=0.5d0*(vec(3,m)+vec(6,m)) 153 vec(4,m)=vec(4,m)-vec(1,m) 154 vec(5,m)=vec(5,m)-vec(2,m) 155 vec(6,m)=vec(6,m)-vec(3,m) 156 endif 157 34 continue 158c 159 do 3 m=1,nmol 160 if(ndx(m,1).gt.0) then 161c 162 tr(1)=-vec(1,m) 163 tr(2)=-vec(2,m) 164 tr(3)=-vec(3,m) 165 if(mode.eq.3) tr(3)=0.0d0 166c 167 do 41 i=1,nmol 168 ndx(i,3)=1 169 41 continue 170 ndx(m,3)=0 171 b2=vec(1,m)*vec(1,m)+vec(2,m)*vec(2,m)+vec(3,m)*vec(3,m) 172 d=sqrt(vec(4,m)*vec(4,m)+vec(5,m)*vec(5,m)+vec(6,m)*vec(6,m)) 173 e=sqrt(vec(4,m)*vec(4,m)+vec(5,m)*vec(5,m)) 174 do 42 i=1,nmol 175 k=1 176 if(ndx(i,1).eq.0.or.i.eq.m) then 177 k=0 178 else 179 if(abs(vec(1,i)-(vec(1,m)+0.5d0*tr(1))).gt. 180 + abs(vec(1,m)+0.5d0*tr(1))+vec(4,i)+vec(4,m)) k=0 181 if(abs(vec(2,i)-(vec(2,m)+0.5d0*tr(2))).gt. 182 + abs(vec(2,m)+0.5d0*tr(2))+vec(5,i)+vec(5,m)) k=0 183 if(abs(vec(3,i)-(vec(3,m)+0.5d0*tr(3))).gt. 184 + abs(vec(3,m)+0.5d0*tr(3))+vec(6,i)+vec(6,m)) k=0 185 if(k.eq.1) then 186 a2=vec(1,i)*vec(1,i)+vec(2,i)*vec(2,i)+vec(3,i)*vec(3,i) 187 c2=(vec(1,m)-vec(1,i))**2+(vec(2,m)-vec(2,i))**2+ 188 + (vec(3,m)-vec(3,i))**2 189 factor=(1.0d0-0.5d0*(b2+c2-a2)/b2) 190 p(1)=factor*vec(1,m) 191 p(2)=factor*vec(2,m) 192 p(3)=factor*vec(3,m) 193 a=sqrt((p(1)-vec(1,i))**2+(p(2)-vec(2,i))**2+(p(3)-vec(3,i))**2) 194 c=sqrt(vec(4,i)**2+vec(5,i)**2+vec(6,i)**2) 195 if(a.gt.c+d) k=0 196 if(mode.eq.3.and. k.eq.1) then 197 if(abs(p(3)-vec(3,i)).gt.abs(vec(6,i))+abs(vec(6,m))) k=0 198 endif 199 endif 200 endif 201 if(k.eq.0) ndx(i,3)=0 202 42 continue 203c 204c determine translation vector 205c 206 rt=-1.0d0 207c 208 dt=sqrt(tr(1)*tr(1)+tr(2)*tr(2)+tr(3)*tr(3)) 209c 210 do 8 i=1,nsa 211 if(ndx(i,2).eq.m) then 212 do 9 j=1,nsa 213 if(ndx(ndx(j,2),3).eq.1) then 214 a=0.0d0 215 b=0.0d0 216 c=0.0d0 217 do 10 k=1,3 218 a=a+tr(k)*tr(k) 219 b=b+(xs(k,i)-xs(k,j))*tr(k) 220 c=c+xs(k,i)*xs(k,i)+xs(k,j)*xs(k,j)-2.0*xs(k,i)*xs(k,j) 221 10 continue 222 b=2.0d0*b 223 c=c-touch*touch 224 if(abs(a).gt.1.0d-5) then 225c 226c find smallest positive root 227c 228 d=(b*b-4.0d0*a*c)/(4.0d0*a*a) 229 if(d.gt.0.0d0) then 230 r1=-b/(2.0d0*a) 231 if(r1.ge.0) then 232 r1=r1+sqrt(d) 233 else 234 r1=r1-sqrt(d) 235 endif 236 r2=c/(r1*a) 237 r=-1.0d0 238 if(r1.gt.0.0d0) r=r1 239 if(r2.gt.0.0d0.and.r2.lt.r1) r=r2 240 if(r.gt.1.0d0) r=-1.0d0 241c 242 if(r.gt.0.0d0) then 243 d=r*sqrt(tr(1)*tr(1)+tr(2)*tr(2)+tr(3)*tr(3)) 244 if(d.lt.dt) then 245 dt=d 246 rt=r 247 endif 248 endif 249c 250 endif 251 endif 252 endif 253 9 continue 254 endif 255 8 continue 256c 257 if(rt.gt.0.and.dt.gt.dtran) then 258 large=m 259 root=rt 260 tran(1)=tr(1) 261 tran(2)=tr(2) 262 tran(3)=tr(3) 263 dtran=dt 264 endif 265 endif 266c 267 3 continue 268c 269 write(*,'(a,i4,a,f12.6,a,f12.6)') 270 + ' Translating molecule ',large,' by ',dtran,' nm',root 271c 272c process largest possible translation 273c 274 if(large.gt.0) then 275 do 11 i=1,nsa 276 if(ndx(i,2).eq.large) then 277 do 12 k=1,3 278 xs(k,i)=xs(k,i)+root*tran(k) 279 12 continue 280 endif 281 11 continue 282 nummov=nummov+1 283 write(filnam,'(a,i5.5,a)') 284 + sysnam(1:index(sysnam,' ')-1)//'-c',nummov,'.pdb' 285c 286 if(.not.pre_wrtpdb(lfnout,lfnpdb,lrgpdb,filnam,iopt,box,num,amass, 287 + mat,csa,isat,isgm,imol,ifra,xs,vs,msa,nsa,cwa,iwat,xw,vw,mwm,mwa, 288 + nwm,nwa,xwc,vwc,mwmc,nwmc,slvnam,iropt,irrand,nxrep,nyrep,nzrep, 289 + drep,msb,nsb,idsb,rdist,nskip,iskip,lang,lfnmrg,nmerge,xmerge, 290 + filmrg,irenum,invert,ihop,ips)) 291 + call md_abort('pre_wrtpdb failed',9999) 292c 293 if(dtran.gt.touch.and.(nummov.lt.nmoves.or.nmoves.lt.0)) goto 2 294 endif 295c 296 return 297 end 298 299c $Id$ 300