1      logical function argos_prep_misfit(lfnout,lfnmat,
2     + isfnd,xs,csa,qsa,isgm,msa,nsa,
3     + idsb,cdsb,msb,nsb,grid,mgrid,ngrid,gdist,nmis,iwater,iwfnd,xw,
4     + qwa,mwm,mwa,nwm,nwa,npbtyp,box,fcount,nrgrid,iogrid,rogrid)
5c
6c $Id$
7c
8      implicit none
9c
10#include "util.fh"
11c
12      logical md_zmat
13      external md_zmat
14c
15      integer lfnout,lfnmat
16      integer msa,nsa,msb,nsb,mgrid,nmis,iwater,mwm,mwa,nwm,nwa
17      integer isfnd(msa),idsb(2,msb),iwfnd(mwm),npbtyp,isgm(msa)
18      real*8 cdsb(msb),xs(3,msa),qsa(msa)
19      character*16 csa(msa)
20      real*8 xw(3,mwa,mwm),qwa(mwa),box(3),fcount
21      character*255 filmat
22      integer nrgrid,iogrid(5),rogrid(2,5)
23c
24      real*8 grid(4,mgrid),distx,disty,distz
25      integer i,ic,j,nn,nndx(10),nnn,nnndx(10),nm,miss
26      real*8 d(3),c(3),p(3),t(3),r,dr,dd,px(3),pn(3),pz(3),pw(3)
27      real*8 xmax(3),xmin(3),dx(3),dist2,gdist,gdist2,angle,boxh(3)
28      integer ix,iy,iz,k,l,m,imax,imin,ngrid,iwm,isa,iwm2,iwa,num
29c
30      integer isg,ioff,jz,kz,lz
31      real*8 dz,az,tz,dgrid
32      logical lmat,fndone
33c
34      if(nmis.gt.0) then
35c
36c     process z-matrices for missing atoms
37c     ------------------------------------
38c
39   32 continue
40      dr=0.0d0
41      ioff=0
42      isg=0
43      lmat=.false.
44      fndone=.false.
45      do 30 i=1,nsa
46      if(isgm(i).ne.isg) then
47      isg=isgm(i)
48      ioff=i-1
49      filmat=csa(i)(1:index(csa(i),' ')-1)//'.mat '
50      lmat=.false.
51      if(lmat) close(unit=lfnmat)
52      open(unit=lfnmat,file=filmat(1:index(filmat,' ')-1),
53     + form='formatted',status='old',err=30)
54      lmat=.true.
55      endif
56      if(lmat.and.isfnd(i).ne.1) then
57      rewind(lfnmat)
58   31 continue
59      read(lfnmat,1000,end=30,err=30) iz,jz,kz,lz,dz,az,tz
60 1000 format(4i5,3f12.6)
61      if(i.eq.iz+ioff.and.isfnd(ioff+jz).ne.0.and.
62     + isfnd(ioff+kz).ne.0.and.isfnd(ioff+lz).ne.0) then
63      if(md_zmat(xs(1,i),xs(1,ioff+jz),xs(1,ioff+kz),xs(1,ioff+lz),
64     + dz,az,tz)) then
65      isfnd(i)=1
66      write(*,'(a,i7,a,a,3f10.5)') ' Generated from z-matrix: ',
67     + isg,':',csa(i),(xs(k,i),k=1,3)
68      fndone=.true.
69      endif
70      else
71      goto 31
72      endif
73      endif
74   30 continue
75      if(lmat) close(unit=lfnmat)
76      if(fndone) goto 32
77c
78      if(util_print('restart',print_default)) then
79      write(lfnout,2000)
80 2000 format(' Z-matrix definitions done',/)
81      endif
82c
83c     generate missing coordinates based on simple rules
84c     --------------------------------------------------
85c
86      gdist2=gdist*gdist
87      boxh(1)=0.5d0*box(1)
88      boxh(2)=0.5d0*box(2)
89      boxh(3)=0.5d0*box(3)
90c
91c     loop over all solute atoms
92c     --------------------------
93c
94      ngrid=0
95      miss=0
96c
97      do 1 i=1,nsa
98c
99c     for each atom without coordinates
100c     ---------------------------------
101c
102      if(isfnd(i).ne.1) then
103      if(util_print('coordinates',print_debug)) then
104      write(lfnout,'(a,i7,5x,a16)') 'Not found: atom ',i,csa(i)
105      endif
106      nn=0
107c
108c     find all bonded atoms
109c     ---------------------
110c
111      do 2 j=1,nsb
112      if(idsb(1,j).eq.i) then
113      nn=nn+1
114      nndx(nn)=idsb(2,j)
115      endif
116      if(idsb(2,j).eq.i) then
117      nn=nn+1
118      nndx(nn)=idsb(1,j)
119      endif
120    2 continue
121c
122      if(nn.eq.0) goto 1
123c
124c     if bound to single other existing atom
125c
126      if(nn.eq.1.and.isfnd(nndx(1)).eq.1) then
127c
128c     find atoms bound to that atom
129c
130      nnn=0
131      do 3 j=1,nsb
132      if(idsb(1,j).eq.nndx(1)) then
133      nnn=nnn+1
134      nnndx(nnn)=idsb(2,j)
135      endif
136      if(idsb(2,j).eq.nndx(1)) then
137      nnn=nnn+1
138      nnndx(nnn)=idsb(1,j)
139      endif
140      if(idsb(1,j).eq.i.or.idsb(2,j).eq.i) dr=cdsb(j)
141    3 continue
142c
143c     count the number of missing atoms
144c
145      nm=0
146      do 4 j=1,nnn
147      if(isfnd(nnndx(j)).ne.1) nm=nm+1
148    4 continue
149c
150c     if j is only missing atom
151c
152      if(nm.eq.1) then
153      c(1)=xs(1,nndx(1))
154      c(2)=xs(2,nndx(1))
155      c(3)=xs(3,nndx(1))
156      d(1)=0.0d0
157      d(2)=0.0d0
158      d(3)=0.0d0
159      num=0
160      do 5 j=1,nnn
161      if(nnndx(j).ne.i) then
162      d(1)=d(1)+xs(1,nnndx(j))
163      d(2)=d(2)+xs(2,nnndx(j))
164      d(3)=d(3)+xs(3,nnndx(j))
165      num=num+1
166      endif
167    5 continue
168      d(1)=d(1)/dble(num)-c(1)
169      d(2)=d(2)/dble(num)-c(2)
170      d(3)=d(3)/dble(num)-c(3)
171      r=sqrt(d(1)*d(1)+d(2)*d(2)+d(3)*d(3))
172      if(r.lt.0.0001)
173     + call md_abort('argos_prep_misfit: vector too small',9001)
174      xs(1,i)=c(1)-dr*d(1)/r
175      xs(2,i)=c(2)-dr*d(2)/r
176      xs(3,i)=c(3)-dr*d(3)/r
177      if(nnn.eq.2) then
178      p(1)=xs(1,i)
179      p(2)=xs(2,i)
180      p(3)=xs(3,i)
181      t(1)=p(2)*c(3)-c(2)*p(3)+c(1)
182      t(2)=p(3)*c(1)-c(3)*p(1)+c(2)
183      t(3)=p(1)*c(2)-c(1)*p(2)+c(3)
184      angle=1.047198d0
185      call rotate(c,t,angle,p,xs(1,i))
186      endif
187      isfnd(i)=1
188      endif
189c
190c     if j is one of two missing atoms
191c
192      if(nm.eq.2) then
193      c(1)=xs(1,nndx(1))
194      c(2)=xs(2,nndx(1))
195      c(3)=xs(3,nndx(1))
196      d(1)=0.0d0
197      d(2)=0.0d0
198      d(3)=0.0d0
199      ic=0
200      num=0
201      do 6 j=1,nnn
202      if(nnndx(j).ne.i.and.isfnd(nnndx(j)).eq.1) then
203      d(1)=d(1)+xs(1,nnndx(j))
204      d(2)=d(2)+xs(2,nnndx(j))
205      d(3)=d(3)+xs(3,nnndx(j))
206      ic=nnndx(j)
207      num=num+1
208      endif
209    6 continue
210      d(1)=d(1)/dble(num)-c(1)
211      d(2)=d(2)/dble(num)-c(2)
212      d(3)=d(3)/dble(num)-c(3)
213      r=sqrt(d(1)*d(1)+d(2)*d(2)+d(3)*d(3))
214      if(r.lt.0.0001)
215     + call md_abort('argos_prep_misfit: vector too small',9002)
216      d(1)=dr*d(1)/r
217      d(2)=dr*d(2)/r
218      d(3)=dr*d(3)/r
219      if(nnn.eq.3) then
220      if(abs(d(1)).le.abs(d(2)).and.abs(d(1)).le.abs(d(3))) then
221      t(1)=d(1)+c(1)
222      t(2)=d(3)+c(2)
223      t(3)=-d(2)+c(3)
224      elseif(abs(d(2)).le.abs(d(1)).and.abs(d(2)).le.abs(d(3))) then
225      t(1)=d(3)+c(1)
226      t(2)=d(2)+c(2)
227      t(3)=-d(1)+c(3)
228      else
229      t(1)=d(2)+c(1)
230      t(2)=-d(1)+c(2)
231      t(3)=d(3)+c(3)
232      endif
233      d(1)=d(1)+c(1)
234      d(2)=d(2)+c(2)
235      d(3)=d(3)+c(3)
236      angle=2.094395102d0
237      call rotate(c,t,angle,d,p)
238      xs(1,i)=p(1)
239      xs(2,i)=p(2)
240      xs(3,i)=p(3)
241      isfnd(i)=1
242      else
243      d(1)=d(1)+c(1)
244      d(2)=d(2)+c(2)
245      d(3)=d(3)+c(3)
246      r=sqrt((xs(1,ic)-c(1))**2+(xs(2,ic)-c(2))**2+(xs(3,ic)-c(3))**2)
247      if(r.lt.0.0001)
248     + call md_abort('argos_prep_misfit: vector too small',9004)
249      t(1)=dr*(xs(1,ic)-c(1))/r+c(1)
250      t(2)=dr*(xs(2,ic)-c(2))/r+c(2)
251      t(3)=dr*(xs(3,ic)-c(3))/r+c(3)
252      angle=1.570796327d0
253      call rotate(c,d,angle,t,p)
254      d(1)=d(1)-c(1)
255      d(2)=d(2)-c(2)
256      d(3)=d(3)-c(3)
257      t(1)=t(1)-c(1)
258      t(2)=t(2)-c(2)
259      t(3)=t(3)-c(3)
260      c(1)=d(2)*t(3)-t(2)*d(3)+xs(1,nndx(1))
261      c(2)=d(3)*t(1)-t(3)*d(1)+xs(2,nndx(1))
262      c(3)=d(1)*t(2)-t(1)*d(2)+xs(3,nndx(1))
263      d(1)=xs(1,nndx(1))
264      d(2)=xs(2,nndx(1))
265      d(3)=xs(3,nndx(1))
266      angle=3.141592654d0
267      call rotate(d,c,angle,p,t)
268      xs(1,i)=t(1)
269      xs(2,i)=t(2)
270      xs(3,i)=t(3)
271      isfnd(i)=1
272      endif
273      endif
274c
275c     if j is one of three missing atoms
276c
277      if(nm.eq.3) then
278      c(1)=xs(1,nndx(1))
279      c(2)=xs(2,nndx(1))
280      c(3)=xs(3,nndx(1))
281      d(1)=0.0d0
282      d(2)=0.0d0
283      d(3)=0.0d0
284      ic=0
285      num=0
286      do 7 j=1,nnn
287      if(nnndx(j).ne.i.and.isfnd(nnndx(j)).eq.1) then
288      d(1)=d(1)+xs(1,nnndx(j))
289      d(2)=d(2)+xs(2,nnndx(j))
290      d(3)=d(3)+xs(3,nnndx(j))
291      ic=nnndx(j)
292      num=num+1
293      endif
294    7 continue
295      d(1)=d(1)/dble(num)-c(1)
296      d(2)=d(2)/dble(num)-c(2)
297      d(3)=d(3)/dble(num)-c(3)
298      r=sqrt(d(1)*d(1)+d(2)*d(2)+d(3)*d(3))
299      if(r.lt.0.0001)
300     + call md_abort('argos_prep_misfit: vector too small',9005)
301      d(1)=dr*d(1)/r
302      d(2)=dr*d(2)/r
303      d(3)=dr*d(3)/r
304      if(nnn.eq.4) then
305      if(abs(d(1)).le.abs(d(2)).and.abs(d(1)).le.abs(d(3))) then
306      t(1)=d(1)+c(1)
307      t(2)=d(3)+c(2)
308      t(3)=-d(2)+c(3)
309      elseif(abs(d(2)).le.abs(d(1)).and.abs(d(2)).le.abs(d(3))) then
310      t(1)=d(3)+c(1)
311      t(2)=d(2)+c(2)
312      t(3)=-d(1)+c(3)
313      else
314      t(1)=d(2)+c(1)
315      t(2)=-d(1)+c(2)
316      t(3)=d(3)+c(3)
317      endif
318      d(1)=d(1)+c(1)
319      d(2)=d(2)+c(2)
320      d(3)=d(3)+c(3)
321      angle=1.910633236d0
322      call rotate(c,t,angle,d,p)
323      xs(1,i)=p(1)
324      xs(2,i)=p(2)
325      xs(3,i)=p(3)
326      isfnd(i)=1
327      else
328c
329c     three missing hydrogens from atom with other than 4 neighbors
330c
331      call md_abort('argos_prep_misfit: not implemented for atom',i)
332      endif
333      endif
334c
335      endif
336      endif
337c
338    1 continue
339c
340      if(util_print('restart',print_default)) then
341      write(lfnout,2001)
342 2001 format(' Geometric rules definitions done',/)
343      endif
344c
345      do 8 i=1,nsa
346c
347c     for each atom without coordinates
348c     ---------------------------------
349c
350      if(isfnd(i).ne.1) then
351      nn=0
352c
353c     find all bonded atoms
354c     ---------------------
355c
356      do 9 j=1,nsb
357      if(idsb(1,j).eq.i) then
358      nn=nn+1
359      nndx(nn)=idsb(2,j)
360      endif
361      if(idsb(2,j).eq.i) then
362      nn=nn+1
363      nndx(nn)=idsb(1,j)
364      endif
365    9 continue
366c
367c
368c     if not bound to any atom
369c     ------------------------
370c
371      if(nn.eq.0) then
372      if(util_print('coordinates',print_debug)) then
373      write(lfnout,'(a,i7)') 'Setting up grid for atom ',i
374      endif
375c
376c     setup grid
377c
378      if(ngrid.eq.0) then
379      do 10 j=1,nsa
380      if(isfnd(j).eq.1) then
381      xmax(1)=xs(1,j)
382      xmax(2)=xs(2,j)
383      xmax(3)=xs(3,j)
384      xmin(1)=xs(1,j)
385      xmin(2)=xs(2,j)
386      xmin(3)=xs(3,j)
387      goto 11
388      endif
389   10 continue
390   11 continue
391      do 12 k=1,3
392      do 13 j=1,nsa
393      if(isfnd(j).eq.1) then
394      if(xmax(k).lt.xs(k,j)) xmax(k)=xs(k,j)
395      if(xmin(k).gt.xs(k,j)) xmin(k)=xs(k,j)
396      endif
397   13 continue
398      dx(k)=(xmax(k)-xmin(k)+4.0d0*gdist)/(mgrid-1)
399      xmin(k)=xmin(k)-2.0d0*gdist
400   12 continue
401c
402      if(ngrid.gt.0.and.util_print('coordinates',print_debug)) then
403      do 555 k=1,nrgrid
404      write(lfnout,'(a,i5,2f12.5)') 'grid option ',
405     + iogrid(k),rogrid(1,k),rogrid(2,k)
406  555 continue
407      endif
408      ngrid=1
409      do 14 ix=1,mgrid
410      do 15 iy=1,mgrid
411      do 16 iz=1,mgrid
412      grid(1,ngrid)=xmin(1)+dble(ix-1)*dx(1)
413      grid(2,ngrid)=xmin(2)+dble(iy-1)*dx(2)
414      grid(3,ngrid)=xmin(3)+dble(iz-1)*dx(3)
415      if(nrgrid.gt.0) then
416      do 217 k=1,nrgrid
417      if(iogrid(k).eq.1) then
418      if(grid(1,ngrid).lt.rogrid(1,k)) goto 16
419      if(grid(1,ngrid).gt.rogrid(2,k)) goto 16
420      elseif(iogrid(k).eq.2) then
421      if(grid(1,ngrid).gt.rogrid(1,k).and.grid(1,ngrid).lt.rogrid(2,k))
422     + goto 16
423      elseif(iogrid(k).eq.3) then
424      if(grid(2,ngrid).lt.rogrid(1,k)) goto 16
425      if(grid(2,ngrid).gt.rogrid(2,k)) goto 16
426      elseif(iogrid(k).eq.4) then
427      if(grid(2,ngrid).gt.rogrid(1,k).and.grid(2,ngrid).lt.rogrid(2,k))
428     + goto 16
429      elseif(iogrid(k).eq.5) then
430      if(grid(3,ngrid).lt.rogrid(1,k)) goto 16
431      if(grid(3,ngrid).gt.rogrid(2,k)) goto 16
432      elseif(iogrid(k).eq.6) then
433      if(grid(3,ngrid).gt.rogrid(1,k).and.grid(3,ngrid).lt.rogrid(2,k))
434     + goto 16
435      elseif(iogrid(k).eq.7) then
436      dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2)
437      if(dgrid.lt.rogrid(1,k)) goto 16
438      if(dgrid.gt.rogrid(2,k)) goto 16
439      elseif(iogrid(k).eq.8) then
440      dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2)
441      if(dgrid.gt.rogrid(1,k).and.dgrid.lt.rogrid(2,k)) goto 16
442      elseif(iogrid(k).eq.9) then
443      dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2+grid(3,ngrid)**2)
444      if(dgrid.lt.rogrid(1,k)) goto 16
445      if(dgrid.gt.rogrid(2,k)) goto 16
446      elseif(iogrid(k).eq.10) then
447      dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2+grid(2,ngrid)**2)
448      if(dgrid.gt.rogrid(1,k).and.dgrid.lt.rogrid(2,k)) goto 16
449      endif
450  217 continue
451      endif
452c      write(*,'(a,4i5)') 'grid ',ix,iy,iz
453      do 17 k=1,nsa
454      if(isfnd(k).eq.1) then
455      if(npbtyp.eq.0) then
456      dist2=(grid(1,ngrid)-xs(1,k))**2+(grid(2,ngrid)-xs(2,k))**2+
457     + (grid(3,ngrid)-xs(3,k))**2
458      else
459      distx=grid(1,ngrid)-xs(1,k)
460      disty=grid(2,ngrid)-xs(2,k)
461      distz=grid(3,ngrid)-xs(3,k)
462      if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1)
463      if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2)
464      if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3)
465      dist2=distx**2+disty**2+distz**2
466      endif
467      if(dist2.lt.gdist2) goto 16
468      endif
469   17 continue
470      if(iwater.gt.0) then
471      do 117 k=1,nwm
472      if(iwfnd(k).gt.0) then
473      do 118 l=1,iwfnd(k)
474      if(npbtyp.eq.0) then
475      dist2=(grid(1,ngrid)-xw(1,l,k))**2+(grid(2,ngrid)-xw(2,l,k))**2+
476     + (grid(3,ngrid)-xw(3,l,k))**2
477      else
478      distx=grid(1,ngrid)-xw(1,l,k)
479      disty=grid(2,ngrid)-xw(2,l,k)
480      distz=grid(3,ngrid)-xw(3,l,k)
481      if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1)
482      if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2)
483      if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3)
484      dist2=distx**2+disty**2+distz**2
485      endif
486      if(dist2.lt.gdist2) goto 16
487  118 continue
488      endif
489  117 continue
490      endif
491      ngrid=ngrid+1
492   16 continue
493   15 continue
494   14 continue
495      ngrid=ngrid-1
496c
497      if(ngrid.le.0) call md_abort('No grid points 0',9999)
498c
499c     calculate grid potential
500c
501      do 18 k=1,ngrid
502      grid(4,k)=0.0d0
503      do 19 j=1,nsa
504      if(isfnd(j).eq.1) then
505      if(npbtyp.eq.0) then
506      dist2=dsqrt((grid(1,k)-xs(1,j))**2+
507     + (grid(2,k)-xs(2,j))**2+(grid(3,k)-xs(3,j))**2)
508      else
509      distx=grid(1,ngrid)-xs(1,k)
510      disty=grid(2,ngrid)-xs(2,k)
511      distz=grid(3,ngrid)-xs(3,k)
512      if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1)
513      if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2)
514      if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3)
515      dist2=distx**2+disty**2+distz**2
516      endif
517      grid(4,k)=grid(4,k)+fcount*qsa(j)/dist2
518      endif
519   19 continue
520      if(iwater.gt.0) then
521      do 119 m=1,nwm
522      if(iwfnd(m).gt.0) then
523      do 120 l=1,iwfnd(m)
524      if(npbtyp.eq.0) then
525      dist2=(grid(1,ngrid)-xw(1,l,m))**2+(grid(2,ngrid)-xw(2,l,m))**2+
526     + (grid(3,ngrid)-xw(3,l,m))**2
527      else
528      distx=grid(1,ngrid)-xw(1,l,m)
529      disty=grid(2,ngrid)-xw(2,l,m)
530      distz=grid(3,ngrid)-xw(3,l,m)
531      if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1)
532      if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2)
533      if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3)
534      dist2=distx**2+disty**2+distz**2
535      endif
536      grid(4,k)=grid(4,k)+qwa(l)/dist2
537  120 continue
538      endif
539  119 continue
540      endif
541   18 continue
542c
543      endif
544c
545      if(ngrid.le.0) call md_abort('No grid points 1',9999)
546c
547c     find extrema
548c
549      imax=1
550      imin=1
551      do 20 j=1,ngrid
552      if(grid(4,j).gt.grid(4,imax)) imax=j
553      if(grid(4,j).lt.grid(4,imin)) imin=j
554   20 continue
555c
556c     choose best grid point
557c
558      if(qsa(i).gt.0.0d0) then
559      xs(1,i)=grid(1,imin)
560      xs(2,i)=grid(2,imin)
561      xs(3,i)=grid(3,imin)
562      else
563      xs(1,i)=grid(1,imax)
564      xs(2,i)=grid(2,imax)
565      xs(3,i)=grid(3,imax)
566      imin=imax
567      endif
568      isfnd(i)=1
569      if(util_print('restart',print_medium)) then
570      write(lfnout,6000) fcount*qsa(i)*grid(4,imin)*138.9354d0
571 6000 format(' Added counterion with energy ',t40,f12.3,' kJ/mol')
572      endif
573c
574c     remove all grid points within gdist from imin from grid
575c
576      d(1)=grid(1,imin)
577      d(2)=grid(2,imin)
578      d(3)=grid(3,imin)
579c
580      imax=ngrid
581      ngrid=0
582      do 21 j=1,imax
583      if(npbtyp.eq.0) then
584      dist2=(grid(1,j)-d(1))**2+(grid(2,j)-d(2))**2+(grid(3,j)-d(3))**2
585      else
586      distx=grid(1,j)-d(1)
587      disty=grid(2,j)-d(2)
588      distz=grid(3,j)-d(3)
589      if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1)
590      if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2)
591      if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3)
592      dist2=distx**2+disty**2+distz**2
593      endif
594      if(dist2.ge.gdist2) then
595      ngrid=ngrid+1
596      grid(1,ngrid)=grid(1,j)
597      grid(2,ngrid)=grid(2,j)
598      grid(3,ngrid)=grid(3,j)
599      dist2=dsqrt((grid(1,ngrid)-xs(1,i))**2+
600     + (grid(2,ngrid)-xs(2,i))**2+(grid(3,ngrid)-xs(3,i))**2)
601      grid(4,ngrid)=grid(4,j)+fcount*qsa(i)/dist2
602      endif
603   21 continue
604c
605      if(ngrid.le.0) call md_abort('No grid points left',9999)
606c
607      else
608      miss=miss+1
609      if(util_print('restart',print_none)) then
610      write(lfnout,9998) isgm(i),csa(i)
611 9998 format(' Coordinates missing for atom ',i5,':',a)
612      endif
613      endif
614      endif
615    8 continue
616c
617      if(util_print('restart',print_default)) then
618      write(lfnout,2002)
619 2002 format(' Grid specifications for non-bonded atoms done',/)
620      endif
621c
622      if(miss.ne.0) then
623      call md_abort('argos_prep_misfit: missing non hydrogen',9999)
624      endif
625c
626      endif
627c
628c     crystal waters
629c
630      if(iwater.gt.0) then
631      do 23 iwm=1,nwm
632      if(iwfnd(iwm).eq.1) then
633      px(1)=0.0d0
634      px(2)=0.0d0
635      px(3)=0.0d0
636      do 24 isa=1,nsa
637      dx(1)=(xs(1,isa)-xw(1,1,iwm))
638      dx(2)=(xs(2,isa)-xw(2,1,iwm))
639      dx(3)=(xs(3,isa)-xw(3,1,iwm))
640      dd=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
641      px(1)=px(1)+qsa(isa)*dx(1)/dd
642      px(2)=px(2)+qsa(isa)*dx(2)/dd
643      px(3)=px(3)+qsa(isa)*dx(3)/dd
644   24 continue
645      do 25 iwm2=1,nwm
646      if(iwfnd(iwm2).eq.3) then
647      do 26 iwa=1,3
648      dx(1)=(xw(1,iwa,iwm2)-xw(1,1,iwm))
649      dx(2)=(xw(2,iwa,iwm2)-xw(2,1,iwm))
650      dx(3)=(xw(3,iwa,iwm2)-xw(3,1,iwm))
651      dd=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
652      px(1)=px(1)+qwa(iwa)*dx(1)/dd
653      px(2)=px(2)+qwa(iwa)*dx(2)/dd
654      px(3)=px(3)+qwa(iwa)*dx(3)/dd
655   26 continue
656      endif
657   25 continue
658      dd=sqrt(px(1)*px(1)+px(2)*px(2)+px(3)*px(3))
659      pw(1)=-0.1d0*px(1)/dd
660      pw(2)=-0.1d0*px(2)/dd
661      pw(3)=-0.1d0*px(3)/dd
662      px(1)=pw(2)
663      px(2)=pw(3)
664      px(3)=pw(1)
665      pn(1)=pw(2)*px(3)-pw(3)*px(2)
666      pn(2)=pw(3)*px(1)-pw(1)*px(3)
667      pn(3)=pw(1)*px(2)-pw(2)*px(1)
668      pz(1)=0.0d0
669      pz(2)=0.0d0
670      pz(3)=0.0d0
671      angle=0.95120d0
672      call rotate(pz,pn,angle,pw,dx)
673      xw(1,2,iwm)=xw(1,1,iwm)+dx(1)
674      xw(2,2,iwm)=xw(2,1,iwm)+dx(2)
675      xw(3,2,iwm)=xw(3,1,iwm)+dx(3)
676      angle=-0.95120d0
677      call rotate(pz,pn,angle,pw,dx)
678      xw(1,3,iwm)=xw(1,1,iwm)+dx(1)
679      xw(2,3,iwm)=xw(2,1,iwm)+dx(2)
680      xw(3,3,iwm)=xw(3,1,iwm)+dx(3)
681c
682      iwfnd(iwm)=3
683      endif
684   23 continue
685      endif
686c
687      if(util_print('restart',print_default)) then
688      write(lfnout,2003)
689 2003 format(' Crystal waters done',/)
690      endif
691c
692      argos_prep_misfit=.true.
693      return
694 9999 continue
695      argos_prep_misfit=.false.
696      return
697      end
698