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