1      subroutine inisw
2      implicit real (a-h,o-z)
3      logical usecut,usesw
4      common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw
5
6      cutof2 = cutoff*cutoff
7      cuton  = cutoff - 1.0e0
8
9      cuton2 = cuton*cuton
10      tmp    = cutof2 - cuton2
11      conof3 = tmp*tmp*tmp
12      isw = 1
13
14
15      return
16      end
17
18      real function swvdw(rij2)
19      implicit real (a-h,o-z)
20      logical usecut,usesw
21      common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw
22
23      if (rij2.ge.cutof2) then
24          swvdw = 0.0e0
25      elseif (rij2.le.cuton2) then
26          swvdw = 1.0e0
27      else
28          rd1 = cutof2 - rij2
29          rd2 = cutof2 + 2.0e0*rij2 - 3.0e0*cuton2
30          swvdw = (rd1*rd1*rd2)/conof3
31      endif
32
33      return
34      end
35
36      real function swdvdw(rij2)
37c
38c derivative of VDW switch function
39c
40      implicit real (a-h,o-z)
41      logical usecut,usesw
42      common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw
43
44      if (rij2.ge.cutof2) then
45          swdvdw = 0.0e0
46      elseif (rij2.le.cuton2) then
47          swdvdw = 1.0e0
48      else
49          rd1 = cutof2 - rij2
50c          rd2 = rij2 - cuton2
51c          swdvdw = 12.0e0*rd1*rd2/conof3
52          rd2 = cutof2 + 2.0e0*rij2 - 3.0e0*cuton2
53          rij = sqrt(rij2)
54          swdvdw = 4.0e0*rij*rd1*(rd1-rd2)/conof3
55      endif
56
57      return
58      end
59
60c swchg not used, instead same function as swvdw used for charge cutoff
61
62      subroutine swchg(rij2,s,ds)
63      implicit real (a-h,o-z)
64      logical usecut,usesw
65      common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw
66
67
68      if (rij2.ge.cutof2) then
69         s = 0.0e0
70         ds = 0.0e0
71      else
72         frac = rij2/cutof2
73         frac2 = frac*frac
74         s = 1.0e0 - frac
75         s = s*s
76         ds = (frac2 - frac)*4.0e0
77c        ds = ds/rij
78c        we put the rij of de = de / r also in here
79c           de = e*ds - e*s/r , de = de / r =>
80c                de = e*ds/rij - e*s/rij2
81c now       de = e*ds - e*s
82         s = s/rij2
83         ds = ds/rij2
84      endif
85
86      return
87      end
88
89      subroutine bldlsd(coo,iresid,nlst,lst,istat)
90      implicit real (a-h,o-z), integer (i-n)
91      parameter (mxneib=200)
92      parameter (numres=50000)
93      logical usecut,usesw
94      common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw
95      common /athlp/  iatoms, mxnat
96      common /residu/  qtot,ihsres,
97     &                 ires(numres),ibeg(numres),iend(numres)
98      logical box,cell,fast
99      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
100      logical doit
101      dimension coo(3,*),iresid(*),nlst(*),lst(*),vr(3),loclst(mxneib)
102
103      if (istat.eq.0) then
104         usesw = .false.
105         print*,'Not enough memory to use cutoff/switch'
106         return
107      endif
108
109      cutf1 = cutoff + 1.0e0
110      cutf2 = cutf1*cutf1
111
112c Build nonbonded interaction list, iatoms*mxneib
113c each atom can have a max of mxneib neighbours within cutoff of 8 angs
114c now not an atom is stored, but a residue
115
116      do i=1,iatoms
117         nlst(i) = 0
118         do j=i+1,iatoms
119
120            if (box) then
121               do k=1,3
122                  vr(k) = coo(k,i) - coo(k,j)
123               end do
124               call reddis(vr)
125               r2 = vr(1)*vr(1) + vr(2)*vr(2) + vr(3)*vr(3)
126               doit = (r2.lt.cutf2)
127            else
128               doit = (dist2(coo(1,i),coo(1,j)).lt.cutf2)
129            endif
130
131            if (doit) then
132               if (nlst(i).lt.mxneib) then
133                  ido = 1
134                  do k=1,nlst(i)
135                     if (iresid(j).eq.lst((i-1)*mxneib+k)) ido = 0
136                  end do
137                  if (ido.eq.1) then
138                     nlst(i) = nlst(i) + 1
139                     lst((i-1)*mxneib+nlst(i)) = iresid(j)
140                  endif
141               else
142                  print*,"neighbour list full for atom ",i
143                  print*,"increase mxneib "
144               endif
145            endif
146         end do
147      end do
148
149      do i=1,ihsres
150         if (ires(i).gt.0) then
151            nloc = nlst(ibeg(i))
152            do j=1,nloc
153               loclst(j) = lst((ibeg(i)-1)*mxneib+j)
154            end do
155            do j=ibeg(i)+1,iend(i)
156               do k=1,nlst(j)
157                  ir = lst((j-1)*mxneib+k)
158                  ido = 1
159                  do l=1,nloc
160                     if (ir.eq.loclst(l)) ido = 0
161                  end do
162                  if (ido.eq.1) then
163                     if (nloc.lt.mxneib) then
164                        nloc = nloc + 1
165                        loclst(nloc) = ir
166                     else
167                        print*,"neighbour list full for atom ",ibeg(i)
168                        print*,"increase mxneib "
169                     endif
170                  endif
171               end do
172            end do
173            do j=ibeg(i),iend(i)
174               nlst(j) = nloc
175               do l=1,nloc
176                  lst((j-1)*mxneib+l) = loclst(l)
177               end do
178            end do
179         endif
180      end do
181
182      return
183      end
184
185      logical function resinc(iatom,nlst,lst,ires)
186      implicit real (a-h,o-z)
187      parameter (mxneib=200)
188      dimension lst(*)
189
190
191      resinc = .false.
192      do i=1,nlst
193         if (lst((iatom-1)*mxneib + i).eq.ires) then
194             resinc = .true.
195             return
196         endif
197      end do
198
199      return
200      end
201
202      real function dist2(a,b)
203c
204c     determine distances between neighboring atoms
205c
206      implicit real (a-h,o-z)
207      dimension a(3)
208      dimension b(3)
209
210      d1 = a(1)-b(1)
211      d2 = a(2)-b(2)
212      d3 = a(3)-b(3)
213
214      dist2 = d1*d1 + d2*d2 + d3*d3
215
216      return
217      end
218
219      subroutine premul
220      implicit real (a-h,o-z), integer (i-n)
221      parameter (mxcls=50)
222      common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls)
223      common /prem/ rs(mxcls,mxcls),es(mxcls,mxcls),iprem
224
225      do i=1,mxcls
226         do j=1,mxcls
227            rst = ambvdwr(i) + ambvdwr(j)
228            rs(i,j) = rst*rst*rst*rst*rst*rst
229            es(i,j) = sqrt(ambvdwe(i)*ambvdwe(j))
230         end do
231      end do
232
233      call allerf
234      call alltab
235
236      iprem = 1
237
238      return
239      end
240
241      subroutine prmulr
242      implicit real(a-h,o-z), integer (i-n)
243      parameter (mxcls=50)
244      real ambvdwr,ambvdwe
245      common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls)
246      common /premr/ rs(mxcls,mxcls),es(mxcls,mxcls),iprem
247
248      do i=1,mxcls
249         do j=1,mxcls
250            rst = real(ambvdwr(i) + ambvdwr(j))
251            rs(i,j) = rst*rst*rst*rst*rst*rst
252            es(i,j) = sqrt(real(ambvdwe(i)*ambvdwe(j)))
253         end do
254      end do
255
256      iprem = 1
257
258      call alltab
259
260      return
261      end
262
263      subroutine alltab
264      implicit real (a-h,o-z), integer (i-n)
265      parameter (mxion=2000)
266      logical box,cell,fast
267      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
268      common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
269      common /atab/ ialtab
270
271      if (ialtab.eq.1) return
272
273      call watcnt
274      call watlst(niwat)
275      call stutab
276
277      ialtab = 1
278
279      return
280      end
281
282      subroutine stutad(woo,whh,woh,dewoo,dewhh,dewoh)
283      implicit real (a-h,o-z), integer (i-n)
284      dimension woo(*),whh(*),woh(*),dewoo(*),dewhh(*),dewoh(*)
285
286      rcut  = 9.0e0
287      rcinv = 1.0e0 / rcut
288      rc2inv =  rcinv * rcinv
289
290      econv = 332.05382e0
291      qo    = -0.8340e0
292      qh    = 0.4170e0
293      vro   = 1.7683e0
294      veo   = 0.1520e0
295      vrh   = 0.0001e0
296      veh   = 0.0e0
297
298      poo   = econv*qo*qo
299      rst   = 2.0e0*vro
300      rst2  = rst*rst
301      pr6oo = rst2*rst2*rst2
302      epsoo = veo
303
304      phh   = econv*qh*qh
305      pr6hh = 0.0e0
306      epshh = 0.0e0
307
308      poh   = econv*qo*qh
309      rst   = vro + vrh
310      rst2  = rst*rst
311      pr6oh = rst2*rst2*rst2
312      epsoh = sqrt(veo*veh)
313
314      nsize = 9000
315
316      do i=1,nsize
317
318         r     = dble(i) / 1000.0e0
319         rinv  = 1.0e0 / r
320         r2inv = rinv*rinv
321         r6inv = r2inv*r2inv*r2inv
322
323         p6    = pr6oo*r6inv
324         p12   = p6*p6
325         eq    = poo * (rinv - rcinv + (r-rcut)*rc2inv)
326         ev    = epsoo*(p12-2.0e0*p6)
327
328         woo(i) = eq + ev
329
330c         if (i.lt.8000) print*,i,' woo ',woo(i),
331c     &                  ' dwoo ',woo(i-1)-woo(i)
332         de = epsoo * (p12 - p6) * (-12.0e0)
333         dewoo(i) =  (de-eq)*r2inv + poo*rinv*rc2inv
334
335         eq    = phh * (rinv - rcinv + (r-rcut)*rc2inv)
336
337         whh(i) = eq
338
339         dewhh(i) =  -eq*r2inv + phh*rinv*rc2inv
340
341         p6    = pr6oh*r6inv
342         p12   = p6*p6
343         eq    = poh * (rinv - rcinv + (r-rcut)*rc2inv)
344         ev    = epsoh*(p12-2.0e0*p6)
345
346         woh(i) = eq + ev
347
348         de = epsoh * (p12 - p6) * (-12.0e0)
349         dewoh(i) =  (de-eq)*r2inv + poh*rinv*rc2inv
350
351      end do
352
353      do i=9000,10000
354         woo(i) = 0.0e0
355         woh(i) = 0.0e0
356         whh(i) = 0.0e0
357         dewoo(i) = 0.0e0
358         dewoh(i) = 0.0e0
359         dewhh(i) = 0.0e0
360      end do
361
362      return
363      end
364
365      subroutine etutad(nsize,woo,whh,woh,dewoo,dewhh,dewoh)
366      implicit real (a-h,o-z), integer (i-n)
367      dimension woo(*),whh(*),woh(*),dewoo(*),dewhh(*),dewoh(*)
368
369      econv = 332.05382e0
370      qo    = -0.8340e0
371      qh    = 0.4170e0
372      vro   = 1.7683e0
373      veo   = 0.1520e0
374      vrh   = 0.0001e0
375      veh   = 0.0e0
376
377      poo   = econv*qo*qo
378      rst   = 2.0e0*vro
379      rst2  = rst*rst
380      pr6oo = rst2*rst2*rst2
381      epsoo = veo
382
383      phh   = econv*qh*qh
384      pr6hh = 0.0e0
385      epshh = 0.0e0
386
387      poh   = econv*qo*qh
388      rst   = vro + vrh
389      rst2  = rst*rst
390      pr6oh = rst2*rst2*rst2
391      epsoh = sqrt(veo*veh)
392
393      do i=1,nsize
394
395         r     = real(i) / 1000.0e0
396         rinv  = 1.0e0 / r
397         r2inv = rinv*rinv
398         r6inv = r2inv*r2inv*r2inv
399
400         p6    = pr6oo*r6inv
401         p12   = p6*p6
402         eq    = poo * rinv
403         ev    = epsoo*(p12-2.0e0*p6)
404
405         woo(i) = eq + ev
406
407         de = epsoo * (p12 - p6) * (-12.0e0)
408         dewoo(i) =  (de-eq)*r2inv
409
410         eq    = phh * rinv
411
412         whh(i) = eq
413
414         dewhh(i) =  -eq*r2inv
415
416         p6    = pr6oh*r6inv
417         p12   = p6*p6
418         eq    = poh * rinv
419         ev    = epsoh*(p12-2.0e0*p6)
420
421         woh(i) = eq + ev
422
423         de = epsoh * (p12 - p6) * (-12.0e0)
424         dewoh(i) =  (de-eq)*r2inv
425
426      end do
427
428      return
429      end
430
431      subroutine allerd(aprerf,aprexp)
432      implicit real (a-h,o-z), integer (i-n)
433      real aprerf,aprexp,alaf,pi,spi,f
434      dimension aprerf(*),aprexp(*)
435
436      alfa   = 0.3e0
437      pi     = 4.e0*atan(1.e0)
438      spi    = sqrt(pi)
439
440      do i=1,1800000
441
442         r     = dble(i) / 1000000.0e0
443         aprerf(i) = real(erfc(r))
444
445      end do
446
447      do i=1,3240000
448
449         f     = - float(i) / 1000000.0e0
450         aprexp(i) = 2.0*alfa*exp(f)/spi
451
452      end do
453
454      return
455      end
456
457c      subroutine allerd(aprerf,aprexp)
458c      implicit real (a-h,o-z), integer (i-n)
459c      real aprerf,aprexp,alaf,pi,spi,f
460c      dimension aprerf(*),aprexp(*)
461c
462c      alfa   = 0.2e0
463c      pi     = 4.e0*atan(1.e0)
464c      spi    = sqrt(pi)
465c
466c      do i=1,3000000
467c
468c         r     = dble(i) / 1000000.0e0
469c         aprerf(i) = real(erfc(r))
470c
471c      end do
472c
473c      do i=1,9000000
474c
475c         f     = - float(i) / 1000000.0e0
476c         aprexp(i) = 2.0*alfa*exp(f)/spi
477c
478c      end do
479c
480c      return
481c      end
482
483      subroutine apperfd(r,fc,aprerf)
484      implicit real (a-h,o-z), integer (i-n)
485      real aprerf,r1000,x,w
486      dimension aprerf(*)
487
488      thou = 1000000.0e0
489
490      if (r.le.10.0e0) then
491          r1000 = real(r)*thou
492          ir = int(r1000)
493          x = (r1000 - float(ir))
494          fc = dble((1.0e0-x)*aprerf(ir) + x*aprerf(ir+1))
495      endif
496
497      return
498      end
499
500