1      subroutine reddis(vr)
2      implicit real (a-h,o-z), integer (i-n)
3      logical box,cell,fast
4      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
5      dimension vr(3)
6
7
8
9      do i=1,3
10         vr(i) = mod(vr(i),abc(i))
11         if (vr(i).gt.abc2(i)) then
12            vr(i) = vr(i) - abc(i)
13         endif
14         if (vr(i).lt.-abc2(i)) then
15            vr(i) = vr(i) + abc(i)
16         endif
17      end do
18
19      return
20      end
21
22      subroutine rddisr(vr)
23      implicit real (a-h,o-z), integer (i-n)
24      real abc,abc2,angles
25      logical box,cell,fast
26      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
27      dimension vr(3)
28
29      do i=1,3
30         do while (vr(i).gt.real(abc2(i)))
31            vr(i) = vr(i) - real(abc(i))
32         end do
33         do while (vr(i).lt.-real(abc2(i)))
34            vr(i) = vr(i) + real(abc(i))
35         end do
36      end do
37
38      return
39      end
40
41      subroutine appbnd(coo,ityp)
42      implicit real (a-h,o-z), integer (i-n)
43      integer*2 ityp
44      logical box,cell,fast,outbox
45      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
46      common /athlp/  iatoms, mxnat
47      dimension coo(3,*),cw(3,3),vec(3),ityp(*)
48
49c put water molecules which are entirely out of the box
50c back in the box
51
52      do i=1,iatoms
53
54         if (ityp(i)  .eq.649.and.
55     &       ityp(i+1).eq.650.and.
56     &       ityp(i+2).eq.650) then
57
58            outbox = .true.
59
60            do j=1,3
61               vec(j) = 0.0e0
62            end do
63
64            do k=1,3
65
66               l = 0
67
68               do j=1,3
69                  cw(j,k) = coo(j,i+k-1)
70                  if (cw(j,k).gt.abc2(j)) then
71                     l = l + 1
72                     if (k.eq.1) vec(j) = -abc(j)
73                  endif
74                  if (cw(j,k).lt.-abc2(j)) then
75                     l = l + 1
76                     if (k.eq.1) vec(j) = abc(j)
77                  endif
78               end do
79
80               if (l.eq.0) outbox = .false.
81
82            end do
83
84c water out of the box, put it back
85
86            if (outbox) then
87               do k=1,3
88                  do j=1,3
89                     coo(j,i+k-1) = cw(j,k) + vec(j)
90                  end do
91               end do
92            endif
93
94         endif
95      end do
96
97      return
98      end
99
100      subroutine makbod(coo)
101      implicit real (a-h,o-z), integer (i-n)
102      logical box,cell,fast
103      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
104      dimension coo(3,*),vec(3)
105
106c     get largest diameter protein and add a default
107c     set abc and abc2
108c     set protein in center of the box
109
110      offs = 14.0e0
111
112      call docnt(vec,coo)
113
114      do i=1,3
115         abc(i) = 2.0e0*vec(i) + offs
116         abc2(i) = 0.5e0*abc(i)
117      end do
118
119      return
120      end
121
122      subroutine allbox
123      implicit real (a-h,o-z), integer (i-n)
124      parameter (mxion=2000)
125      logical box,cell,fast
126      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
127      common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
128      common /athlp/  iatoms, mxnat
129      dimension nbox(3)
130
131      watd = 18.620e0
132      nwater = 648/3
133
134      do i=1,3
135         t = abc(i)/watd
136         nbox(i) = int(t) + 1
137      end do
138
139      newat = nbox(1)*nbox(2)*nbox(3)*nwater*3
140
141      call allcoo(newat)
142      if (nion.gt.0) call allq
143
144      return
145      end
146
147      subroutine filbod(water,coo,iconn,iresid,ityp,
148     &                  nac,iac,nad,
149     &                  nbnd,ibnd,bl,bk,
150     &                  nang,iang,ango,ak,q,iopt,iwtpr)
151      implicit real (a-h,o-z), integer (i-n)
152      parameter (mxcon=10)
153      parameter (mxac=3*mxcon)
154      parameter (numres=50000)
155      parameter (mxion=2000)
156      logical box,cell,fast,chkwat,chkbox,chkion
157      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
158      common /athlp/  iatoms, mxnat
159      common /residu/  qtot,ihsres,
160     &                 ires(numres),ibeg(numres),iend(numres)
161      common /h2oer/   numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
162      integer*2 ityp
163      dimension water(3,*), coo(3,*), nbox(3), doff(3)
164      dimension ityp(*),iconn(mxcon+1,*),iresid(*)
165      dimension nac(*),nad(*),iac(mxac,*)
166      dimension ibnd(2,*), bl(*), bk(*)
167      dimension iang(3,*),ango(*),ak(*)
168      dimension q(*),iopt(*),iwtpr(*),cwat(3,3),ibox(3)
169
170      watd = 18.620e0
171      nwater = 648/3
172
173c 18.620 ang box
174c     take prefilled water box and tile the real box with this
175c     sub box
176c     remove water that overlap with protein
177
178      do i=1,3
179         t = abc(i)/watd
180         nbox(i) = int(t) + 1
181      end do
182
183      ishoh = ires(ihsres)
184      if (ishoh.gt.0) then
185          ishoh = -4
186      else
187          ishoh = ishoh - 1
188      endif
189
190      natnow = iatoms
191      ioff = iatoms
192      numwat = 0
193      nwat = 0
194
195      ibox(1) = 0
196
197      do i=1,nbox(1)
198         if (i.eq.nbox(1)) ibox(1) = 1
199         ibox(2) = 0
200
201         do j=1,nbox(2)
202            if (j.eq.nbox(2)) ibox(2) = 1
203            ibox(3) = 0
204
205            do k=1,nbox(3)
206               if (k.eq.nbox(3)) ibox(3) = 1
207
208               doff(1) = -abc2(1) + 0.5e0*watd + watd*dble(i-1)
209               doff(2) = -abc2(2) + 0.5e0*watd + watd*dble(j-1)
210               doff(3) = -abc2(3) + 0.5e0*watd + watd*dble(k-1)
211
212               chkbox = .false.
213               if (i.eq.nbox(1).or.j.eq.nbox(2).or.k.eq.nbox(3))
214     &             chkbox = .true.
215
216               do n=1,nwater
217
218                  in = (n-1)*3
219
220                  do m=1,3
221                     cwat(m,1) = water(m,in+1) + doff(m)
222                     cwat(m,2) = water(m,in+2) + doff(m)
223                     cwat(m,3) = water(m,in+3) + doff(m)
224                  end do
225
226                  if (chkwat(cwat,coo,ityp,iopt,ioff,ibox,
227     &                       iwt,chkbox)) then
228
229
230                     nwat = nwat + 1
231
232                     if (.not.chkion(nwat)) then
233
234c ADD WATER
235
236                        numwat = numwat + 1
237                        iwtpr(numwat) = iwt
238
239                        if (ihsres.lt.numres) then
240                           ihsres = ihsres + 1
241                        endif
242
243                        ires(ihsres) = ishoh
244                        ibeg(ihsres) = ioff+1
245                        iend(ihsres) = ioff+3
246                        iresid(ioff+1) = ishoh
247                        iresid(ioff+2) = ishoh
248                        iresid(ioff+3) = ishoh
249                        ishoh = ishoh - 1
250
251                        do m=1,3
252                           coo(m,ioff+1) = cwat(m,1)
253                           coo(m,ioff+2) = cwat(m,2)
254                           coo(m,ioff+3) = cwat(m,3)
255                        end do
256
257                        ityp(ioff+1) = 649
258                        ityp(ioff+2) = 650
259                        ityp(ioff+3) = 650
260
261                        iwtpr(numwat) = 0
262
263                        iconn(1,ioff+1) = 2
264                        iconn(2,ioff+1) = ioff+2
265                        iconn(3,ioff+1) = ioff+3
266
267                        iconn(1,ioff+2) = 1
268                        iconn(2,ioff+2) = ioff+1
269
270                        iconn(1,ioff+3) = 1
271                        iconn(2,ioff+3) = ioff+1
272
273                        ibnd(1,nbnd+1) = ioff+1
274                        ibnd(2,nbnd+1) = ioff+2
275                        bl(nbnd+1)     = 0.9572e0
276                        bk(nbnd+1)     = 553.0e0
277
278                        ibnd(1,nbnd+2) = ioff+1
279                        ibnd(2,nbnd+2) = ioff+3
280                        bl(nbnd+2)     = 0.9572e0
281                        bk(nbnd+2)     = 553.0e0
282                        nbnd = nbnd + 2
283
284                        nang = nang + 1
285                        iang(1,nang) = ioff+2
286                        iang(2,nang) = ioff+1
287                        iang(3,nang) = ioff+3
288                        ango(nang)   = 104.52e0
289                        ak(nang)     = 100.00e0
290
291                        q(ioff+1)    = -0.8340e0
292                        q(ioff+2)    = 0.4170e0
293                        q(ioff+3)    = 0.4170e0
294
295                        if (i.eq.1.or.j.eq.1.or.k.eq.1) then
296                           iopt(ioff+1) = 1
297                           iopt(ioff+2) = 1
298                           iopt(ioff+3) = 1
299                        else
300                           iopt(ioff+1) = 0
301                           iopt(ioff+2) = 0
302                           iopt(ioff+3) = 0
303                        endif
304
305                        nac(ioff+1) = 0
306                        nac(ioff+2) = 1
307                        nac(ioff+3) = 1
308                        iac(1,ioff+2) = ioff+3
309                        iac(1,ioff+3) = ioff+2
310                        nad(ioff+1) = 0
311                        nad(ioff+2) = 0
312                        nad(ioff+3) = 0
313
314                        ioff = ioff + 3
315
316                     endif
317
318                  endif
319
320               end do
321
322            end do
323         end do
324      end do
325
326      do i=iatoms,ioff
327         iopt(i) = 1
328      end do
329
330      iatoms = ioff
331
332      if (nion.ne.0)
333     &  call addions(coo,iconn,iresid,ityp,q,iopt)
334
335      return
336      end
337
338      logical function chkion(iwat)
339      implicit real (a-h,o-z), integer (i-n)
340      parameter (mxion=2000)
341      common /h2oer/  numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
342
343      chkion = .false.
344      if (nion.eq.0) return
345
346      do i=1,nion
347         if (iwat.eq.iw(i)) then
348            chkion = .true.
349            return
350         endif
351      end do
352
353      return
354      end
355
356      subroutine addions(coo,iconn,iresid,ityp,q,iopt)
357      implicit real (a-h,o-z), integer (i-n)
358      parameter (mxcon=10)
359      parameter (numres=50000)
360      parameter (mxion=2000)
361      common /athlp/  iatoms, mxnat
362      common /residu/  qtot,ihsres,
363     &                 ires(numres),ibeg(numres),iend(numres)
364      common /h2oer/   numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
365      integer*2 ityp
366      dimension coo(3,*), ityp(*),iconn(mxcon+1,*),iresid(*)
367      dimension q(*),iopt(*)
368
369c ADD ION
370
371      ioff = iatoms
372      ishoh = iresid(ioff) - 1
373
374      do i=1,nion
375         ihsres = ihsres + 1
376         ires(ihsres) = ishoh
377         ibeg(ihsres) = ioff+1
378         iend(ihsres) = ioff+1
379         iresid(ioff+1) = ishoh
380
381         ishoh = ishoh - 1
382
383         iatptr = natnow + (iw(i)-1)*3 + 1
384
385         do j=1,3
386            coo(j,ioff+1) = coo(j,iatptr)
387         end do
388
389         if (iontyp.eq.1) then
390            ityp(ioff+1) = 659
391            q(ioff+1)    = -1.0e0
392         else
393            ityp(ioff+1) = 652
394            q(ioff+1)    = +1.0e0
395         endif
396
397         iconn(1,ioff+1) = 0
398         iopt(ioff+1) = 1
399
400         ioff = ioff + 1
401
402      end do
403
404      iatoms = ioff
405
406      return
407      end
408
409      subroutine cntwad(water,coo,iconn,iresid,ityp,q,iopt,qpot)
410      implicit real (a-h,o-z), integer (i-n)
411      parameter (mxcon=10)
412      parameter (mxac=3*mxcon)
413      parameter (numres=50000)
414      parameter (mxion=2000)
415      logical doind,exclu,box,cell,fast,chkwat,chkbox
416      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
417      common /athlp/  iatoms, mxnat
418      common /residu/  qtot,ihsres,
419     &                 ires(numres),ibeg(numres),iend(numres)
420      common /h2oer/   numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
421      integer*2 ityp
422      dimension water(3,*), coo(3,*), nbox(3), doff(3)
423      dimension ityp(*),iconn(mxcon+1,*),iresid(*)
424      dimension q(*),iopt(*),cwat(3,3),c1(3),c2(3),ibox(3),qpot(*)
425
426      watd = 18.620e0
427      nwater = 648/3
428
429c 18.620 ang box
430c     take prefilled water box and tile the real box with this
431c     sub box
432c     remove water that overlap with protein
433
434      do i=1,3
435         t = abc(i)/watd
436         nbox(i) = int(t) + 1
437      end do
438
439      ioff = iatoms
440      nitmp = 0
441
442      ibox(1) = 0
443
444      do i=1,nbox(1)
445         if (i.eq.nbox(1)) ibox(1) = 1
446         ibox(2) = 0
447
448         do j=1,nbox(2)
449            if (j.eq.nbox(2)) ibox(2) = 1
450            ibox(3) = 0
451
452            do k=1,nbox(3)
453               if (k.eq.nbox(3)) ibox(3) = 1
454
455               doff(1) = -abc2(1) + 0.5e0*watd + watd*dble(i-1)
456               doff(2) = -abc2(2) + 0.5e0*watd + watd*dble(j-1)
457               doff(3) = -abc2(3) + 0.5e0*watd + watd*dble(k-1)
458
459               chkbox = .false.
460               if (i.eq.nbox(1).or.j.eq.nbox(2).or.k.eq.nbox(3))
461     &             chkbox = .true.
462
463               do n=1,nwater
464
465                  in = (n-1)*3
466
467                  do m=1,3
468                     cwat(m,1) = water(m,in+1) + doff(m)
469                     cwat(m,2) = water(m,in+2) + doff(m)
470                     cwat(m,3) = water(m,in+3) + doff(m)
471                  end do
472
473                  if (chkwat(cwat,coo,ityp,iopt,ioff,ibox,
474     &                       iwt,chkbox)) then
475
476                     numwat = numwat + 1
477
478                     if (i.eq.1.or.j.eq.1.or.k.eq.1) then
479                           iopt(ioff+1) = 1
480                           iopt(ioff+2) = 1
481                           iopt(ioff+3) = 1
482                     else
483                           iopt(ioff+1) = 0
484                           iopt(ioff+2) = 0
485                           iopt(ioff+3) = 0
486                     endif
487
488                     ityp(ioff+1) = 649
489                     ityp(ioff+2) = 650
490                     ityp(ioff+3) = 650
491
492                     do m=1,3
493                        coo(m,ioff+1) = cwat(m,1)
494                        coo(m,ioff+2) = cwat(m,2)
495                        coo(m,ioff+3) = cwat(m,3)
496                     end do
497
498                     if (ionpl.eq.1) then
499                        call clmond(cwat(1,1),pot,coo,q)
500                        qpot(numwat) = pot
501                     endif
502
503                     ioff = ioff + 3
504
505                  endif
506
507               end do
508
509            end do
510         end do
511      end do
512
513      if (ionpl.eq.1) then
514
515         nitmp = 0
516
517         do while (nitmp.lt.nion)
518
519            doind = .false.
520            potl = 1.0e10
521            iptr = 0
522
523            do i=1,numwat
524
525               pot = qpot(i)
526
527               exclu = .false.
528               do j=1,nitmp
529                  if (iw(j).eq.i) exclu = .true.
530               end do
531
532               if (pot.lt.potl.and..not.exclu) then
533                  doind = .true.
534                  potl = pot
535                  iptr = i
536               endif
537
538            end do
539
540            if (doind) then
541               nitmp = nitmp + 1
542               iw(nitmp) = iptr
543
544               do i=1,numwat
545                  iatptr = iatoms + (i-1)*3 + 1
546                  do j=1,3
547                     c1(j) = coo(j,iatptr)
548                  end do
549
550                  jatptr = iatoms + (iptr-1)*3 + 1
551                  do j=1,3
552                     c2(j) = coo(j,jatptr)
553                  end do
554                  r2 = dist2(c1,c2)
555                  r = sqrt(r2)
556                  qpot(i) = qpot(i) + 1.0e0/r
557
558               end do
559
560            endif
561
562         end do
563
564      else
565
566         call crionp
567
568      endif
569
570      return
571      end
572
573      logical function chkvec(vec,cwat,coo,iopt,ityp,ioff)
574      implicit real (a-h,p-w),integer (i-n),logical (o)
575      parameter (mxamb=1590)
576      parameter (mxgff=72)
577      integer ambcls
578      common /typpar/ ambchg(mxamb),ambcls(mxamb)
579      parameter (mxcls=49)
580      common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls)
581      integer amb2gf
582      common /gftyp/  gfvdw(2,mxgff),amb2gf(mxamb)
583      common /athlp/  iatoms, mxnat
584      integer*2 ityp
585      dimension v(3),vec(3),cwat(3,3),coo(3,*),ityp(*)
586      dimension iopt(*)
587
588      chkvec = .true.
589
590      vdwat = 2.7683e0
591
592      do i=iatoms,ioff
593
594         if (iopt(i).eq.1) then
595
596            i1 = int(ityp(i))
597            if (i1.gt.0) then
598               il = ambcls(i1)
599               vdwr = ambvdwr(il)
600            else
601               il = iabs(i1)
602               vdwr = gfvdw(1,il)
603            endif
604
605
606               do k=1,3
607                  v(k) = coo(k,i) + vec(k) - cwat(k,1)
608               end do
609
610               rab2 = v(1)*v(1) + v(2)*v(2) + v(3)*v(3)
611               dmaxsq = vdwr + vdwat
612               dmaxsq = dmaxsq * dmaxsq
613
614               if (rab2.lt.dmaxsq) then
615                  chkvec = .false.
616                  return
617               endif
618
619
620         endif
621
622      end do
623
624      return
625      end
626
627      logical function chkwat(cwat,coo,ityp,iopt,ioff,ibox,iwrpr,chkbox)
628      implicit real (a-h,p-w),integer (i-n),logical (o)
629      parameter (mxamb=1590)
630      integer ambcls
631      common /typpar/ ambchg(mxamb),ambcls(mxamb)
632      parameter (mxcls=49)
633      common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls)
634      parameter (mxgff=72)
635      integer amb2gf
636      common /gftyp/  gfvdw(2,mxgff),amb2gf(mxamb)
637      integer*2 ityp
638      logical chkbox,chkvec
639      common /athlp/  iatoms, mxnat
640      logical box,cell,fast,owat
641      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
642      dimension coo(3,*),ityp(*),cwat(3,3),v(3),vdwat(3),ibox(3)
643      dimension vec(3),iopt(*)
644
645      chkwat = .false.
646
647
648      owat = .true.
649      iwrpr = 0
650
651      if (chkbox) then
652
653         do i=1,3
654            do j=1,3
655               if (cwat(j,i).gt.abc2(j)) owat = .false.
656            end do
657         end do
658
659         if (.not.owat) return
660
661         do i=1,3
662
663            do j=1,3
664               vec(j) = 0.0e0
665            end do
666
667            if (ibox(i).eq.1) then
668               vec(i) = abc(i)
669               if (.not.chkvec(vec,cwat,coo,iopt,ityp,ioff)) return
670            endif
671
672         end do
673
674c         chkwat = .true.
675c         return
676
677      endif
678
679c check overlap with protein atoms
680
681      owat = .true.
682
683      vdwat(1) = 2.7683e0
684
685      do i=1,iatoms
686
687         i1 = int(ityp(i))
688         if (i1.gt.0) then
689            il = ambcls(i1)
690            vdwr = ambvdwr(il)
691         elseif (i1.le.0) then
692            i1 = iabs(i1)
693            vdwr = gfvdw(1,i1)
694         endif
695
696         do k=1,3
697            v(k) = coo(k,i) - cwat(k,1)
698         end do
699         rab2 = v(1)*v(1) + v(2)*v(2) + v(3)*v(3)
700         dmaxsq = vdwr + vdwat(1)
701         dmaxsq = dmaxsq * dmaxsq
702         if (rab2.lt.dmaxsq) owat = .false.
703         if (rab2.lt.81.0e0) iwrpr = 1
704
705      end do
706
707      if (owat) chkwat = .true.
708
709      return
710      end
711
712      subroutine docnt(vecm,coo)
713      implicit real (a-h,p-w),integer (i-n),logical (o)
714      common /athlp/ iatoms, mxnat
715      dimension vec(3),vecm(3)
716      dimension coo(3,*)
717
718      do i=1,3
719         vecm(i) = 0.0e0
720      end do
721
722      call cntvec(vec,coo,iatoms)
723
724      do i=1,iatoms
725         do j=1,3
726            coo(j,i) = coo(j,i) - vec(j)
727            da = abs(coo(j,i))
728            if (da.gt.vecm(j)) vecm(j) = da
729         end do
730      end do
731
732      return
733      end
734
735      subroutine cntvec(vec,coo,iatoms)
736      implicit real (a-h,p-w),integer (i-n),logical (o)
737      dimension vec(3), coo(3,*)
738
739      do i=1,3
740         vec(i) = 0.0e0
741      end do
742
743      if (iatoms.le.0) return
744
745      do i=1,iatoms
746         do j=1,3
747            vec(j) = vec(j) + coo(j,i)
748         end do
749      end do
750
751      do j=1,3
752         vec(j) = vec(j) / dble(iatoms)
753      end do
754
755      return
756      end
757
758      subroutine setop(a,b,c,alpha,beta,gamma)
759      implicit real (a-h,p-z),integer (i-n),logical (o)
760      common /cell/ xa,ya,yb,za,zb,zc
761      real rxa,rya,ryb,rza,rzb,rzc
762      common /cellr/ rxa,rya,ryb,rza,rzb,rzc
763
764      torad = atan(1.0e0) / 45.0e0
765      alpha = alpha*torad
766      beta = beta*torad
767      gamma = gamma*torad
768
769      ca = cos(alpha)
770      cb = cos(beta)
771      cc = cos(gamma)
772      sc = sin(gamma)
773      sa = sin(alpha)
774      cvol = a*b*c*
775     &   sqrt(1-ca**2-cb**2-cc**2+2.0e0*ca*cb*cc)
776      xa = cvol / (b*c*sa)
777      ya = (a*(cc-(ca*cb)))/sa
778      yb = (b*sa)
779      za = (a*cb)
780      zb = (b*ca)
781      zc = c
782
783      rxa = real(xa)
784      rya = real(ya)
785      ryb = real(yb)
786      rza = real(za)
787      rzb = real(zb)
788      rzc = real(zc)
789
790c      print*,' '
791c      write(*,'(a,f9.3)') ' Cell Volume = ',cvol
792
793      return
794      end
795
796      subroutine getcel(a,b,c,alpha,beta,gamma)
797      implicit real (a-h,p-z),integer (i-n),logical (o)
798      common /cell/ xa,ya,yb,za,zb,zc
799
800      torad = atan(1.0e0) / 45.0e0
801
802      a = sqrt(xa**2+ya**2+za**2)
803      b = sqrt(yb**2+zb**2)
804      c = zc
805
806      sa = yb/b
807      ca = zb/b
808      cb = za/a
809      cc = ca*cb + ya*sa/a
810
811      alpha = acos(ca) / torad
812      beta  = acos(cb) / torad
813      gamma = acos(cc) / torad
814
815      return
816      end
817
818      subroutine fr2crt(vec)
819      implicit real (a-h,p-z),integer (i-n),logical (o)
820      common /cell/ xa,ya,yb,za,zb,zc
821      dimension vec(3)
822
823      vec(3) = za*vec(1) + zb*vec(2) + zc*vec(3)
824      vec(2) = ya*vec(1) + yb*vec(2)
825      vec(1) = xa*vec(1)
826
827      return
828      end
829
830      subroutine crt2fr(veci,veco)
831      implicit real (a-h,p-z),integer (i-n),logical (o)
832      common /cell/ xa,ya,yb,za,zb,zc
833      dimension veci(3),veco(3)
834
835      veco(1) =  veci(1) / xa
836      veco(2) = (veci(2) - ya*veco(1)) / yb
837      veco(3) = (veci(3) - za*veco(1) - zb*veco(2)) / zc
838
839      return
840      end
841
842      subroutine fr2crr(vec)
843      implicit real (a-h,p-z),integer (i-n),logical (o)
844      common /cellr/ xa,ya,yb,za,zb,zc
845      dimension vec(3)
846
847      vec(3) = za*vec(1) + zb*vec(2) + zc*vec(3)
848      vec(2) = ya*vec(1) + yb*vec(2)
849      vec(1) = xa*vec(1)
850
851      return
852      end
853
854      subroutine crr2fr(veci,veco)
855      implicit real (a-h,p-z),integer (i-n),logical (o)
856      common /cellr/ xa,ya,yb,za,zb,zc
857      dimension veci(3),veco(3)
858
859      veco(1) =  veci(1) / xa
860      veco(2) = (veci(2) - ya*veco(1)) / yb
861      veco(3) = (veci(3) - za*veco(1) - zb*veco(2)) / zc
862
863      return
864      end
865
866      subroutine redfr(xyz)
867      implicit real (a-h,p-z),integer (i-n),logical (o)
868      logical box,cell,fast
869      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
870      dimension xyz(3),fr(3)
871
872      call crt2fr(xyz,fr)
873
874      do i=1,3
875         frt = 1.0e0
876         if (fr(i).lt.0.0e0) then
877             fr(i) = abs(fr(i))
878             frt = -1.0e0
879         endif
880
881         fr(i) = mod(fr(i),1.0e0)
882         if (fr(i) .gt. 0.5e0) fr(i) = fr(i) - 1.0e0
883
884         xyz(i) = frt*fr(i)
885      end do
886
887      call fr2crt(xyz)
888
889      return
890      end
891
892      subroutine rddfrr(xyz)
893      implicit real (a-h,p-z),integer (i-n),logical (o)
894      real abc,abc2,angles
895      logical box,cell,fast
896      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
897      common /debug/ ideb
898      dimension xyz(3),fr(3)
899
900
901      call crr2fr(xyz,fr)
902
903
904      do i=1,3
905         frt = 1.0e0
906         if (fr(i).lt.0.0e0) then
907             fr(i) = abs(fr(i))
908             frt = -1.0e0
909         endif
910
911         fr(i) = amod(fr(i),1.0e0)
912         if (fr(i) .gt. 0.5e0) fr(i) = fr(i) - 1.0e0
913
914         xyz(i) = frt*fr(i)
915      end do
916
917
918      call fr2crr(xyz)
919
920      return
921      end
922
923      subroutine adcldr(cellder,ded,de,fct,lex)
924      implicit real (a-h,p-z),integer (i-n),logical (o)
925      common /cell/ xa,ya,yb,za,zb,zc
926      common /expnd/ addfr(3,125),ie(5),izero
927      dimension ded(3),cellder(6),addf(3)
928
929      do j=1,3
930         addf(j) = addfr(j,lex)
931      end do
932
933      cellder(1) = cellder(1) + addf(1)*ded(1)*de*fct
934
935      cellder(2) = cellder(2) + addf(1)*ded(2)*de*fct
936      cellder(3) = cellder(3) + addf(2)*ded(2)*de*fct
937
938      cellder(4) = cellder(4) + addf(1)*ded(3)*de*fct
939      cellder(5) = cellder(5) + addf(2)*ded(3)*de*fct
940      cellder(6) = cellder(6) + addf(3)*ded(3)*de*fct
941
942      return
943      end
944
945      subroutine prcldr(cellder,par)
946      implicit real (a-h,p-z),integer (i-n),logical (o)
947      common /athlp/ iatoms, mxnat
948      dimension cellder(6),par(3,*)
949
950c      do i=1,6
951c         print*,'cellder (',i,')=',cellder(i)
952c      end do
953
954      print*,'val ',par(1,iatoms+1),' cellder (1)=',cellder(1)
955      print*,'val ',par(2,iatoms+1),' cellder (2)=',cellder(2)
956      print*,'val ',par(3,iatoms+1),' cellder (3)=',cellder(3)
957      print*,'val ',par(1,iatoms+2),' cellder (4)=',cellder(4)
958      print*,'val ',par(2,iatoms+2),' cellder (5)=',cellder(5)
959      print*,'val ',par(3,iatoms+2),' cellder (6)=',cellder(6)
960
961      return
962      end
963
964      subroutine expfr(in,cxyz,exyz,lex,nexpnd,oje)
965      implicit real (a-h,p-z),integer (i-n),logical (o)
966      common /expnd/ addfr(3,125),ie(5),izero
967      logical oje(125)
968      dimension lex(125),cxyz(3,125),exyz(3,125,*)
969
970      l = 0
971
972      do i=1,125
973
974         if (i.ne.izero) then
975            l = l + 1
976            do j=1,3
977               cxyz(j,l) = exyz(j,i,in) - exyz(j,izero,in)
978            end do
979            lex(l) = i
980            oje(l) = .false.
981         endif
982
983      end do
984
985      return
986      end
987
988      subroutine setfrac(coo,frac,exyz)
989      implicit real (a-h,p-z),integer (i-n),logical (o)
990      common /athlp/  iatoms, mxnat
991      common /expnd/ addfr(3,125),ie(5),izero
992      dimension coo(3,*),frac(3,*),exyz(3,125,*),cxyz(3,125),dfr(3)
993
994      do j=1,3
995         dfr(j) = 0.0e0
996      end do
997
998      do i=1,iatoms
999
1000         do j=1,3
1001            if (abs(coo(j,i)).lt.0.001e0) coo(j,i) = 0.0e0
1002         end do
1003
1004         call crt2fr(coo(1,i),frac(1,i))
1005
1006         do j=1,3
1007            if (frac(j,i).lt.dfr(j)) dfr(j) = frac(j,i)
1008         end do
1009
1010      end do
1011
1012c huh ???
1013
1014      do i=1,5
1015         do j=1,3
1016            if ((dfr(j).ge.dble(-1*i)).and.(dfr(j).lt.dble(-1*(i-1))))
1017     &          dfr(j) = dble(-1*i)
1018         end do
1019      end do
1020
1021      do i=1,iatoms
1022
1023          do j=1,3
1024             frac(j,i) = frac(j,i) - dfr(j)
1025          end do
1026
1027          do k=1,125
1028             do j=1,3
1029                exyz(j,k,i) = frac(j,i) + addfr(j,k)
1030             end do
1031             call fr2crt(exyz(1,k,i))
1032          end do
1033
1034      end do
1035
1036      return
1037      end
1038
1039      subroutine xpfr(in,kn,cxyz,exyz,lex,nexpnd,oje)
1040      implicit real (a-h,p-z),integer (i-n),logical (o)
1041      common /expnd/ addfr(3,125),ie(5),izero
1042      logical oje(125)
1043      dimension cxyz(3,125),exyz(3,125,*),lex(125)
1044
1045      do i=1,nexpnd+1
1046         oje(i) = .false.
1047         if (i.eq.izero) oje(i) = .true.
1048         do j=1,3
1049            cxyz(j,i) = exyz(j,i,in) - exyz(j,izero,kn)
1050         end do
1051         lex(i) = i
1052      end do
1053
1054      return
1055      end
1056
1057      subroutine initxp
1058      implicit real (a-h,p-z),integer (i-n),logical (o)
1059      common /expnd/ addfr(3,125),ie(5),izero
1060      data (ie(j),j=1,5)  / -2,-1,0,1,2/
1061
1062      izero = 0
1063
1064      l = 0
1065      do i=1,5
1066         do j=1,5
1067            do k=1,5
1068              l = l + 1
1069              addfr(1,l) = dble(ie(i))
1070              addfr(2,l) = dble(ie(j))
1071              addfr(3,l) = dble(ie(k))
1072              if (ie(i).eq.0.and.ie(j).eq.0.and.ie(k).eq.0) izero = l
1073            end do
1074         end do
1075      end do
1076
1077      return
1078      end
1079
1080      subroutine var2cl(x,n)
1081c copy cell variables to cell common
1082      implicit real (a-h,p-z),integer (i-n),logical (o)
1083      logical box,cell,fast
1084      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
1085      common /cell/ xa,ya,yb,za,zb,zc
1086      real rxa,rya,ryb,rza,rzb,rzc
1087      common /cellr/ rxa,rya,ryb,rza,rzb,rzc
1088      dimension x(*)
1089
1090      n3 = n*3
1091
1092      xa = x(n3+1)
1093      ya = x(n3+2)
1094      yb = x(n3+3)
1095      za = x(n3+4)
1096      zb = x(n3+5)
1097      zc = x(n3+6)
1098
1099      abc(1) = sqrt(xa**2+ya**2+za**2)
1100      abc(2) = sqrt(yb**2+zb**2)
1101      abc(3) = zc
1102
1103      do i=1,3
1104         abc2(i) = abc(i)*0.5e0
1105      end do
1106
1107      rxa = real(x(n3+1))
1108      rya = real(x(n3+2))
1109      ryb = real(x(n3+3))
1110      rza = real(x(n3+4))
1111      rzb = real(x(n3+5))
1112      rzc = real(x(n3+6))
1113
1114      return
1115      end
1116
1117      subroutine grd2var(cellder,f,n)
1118c copy cell variables to cell common
1119      implicit real (a-h,p-z),integer (i-n),logical (o)
1120      common /cell/ xa,ya,yb,za,zb,zc
1121      dimension cellder(6),f(*)
1122
1123c      print*,'cellval ',xa,ya,yb,za,zb,zc
1124
1125      n3 = n*3
1126
1127      do i=1,6
1128         f(n3+i) = -cellder(i)
1129      end do
1130
1131      return
1132      end
1133
1134      subroutine uinner(uin,cellder,coo,f,q,n)
1135      implicit real (a-h,p-z),integer (i-n),logical (o)
1136      common /cell/ xa,ya,yb,za,zb,zc
1137      dimension cellder(6),pole(3),coo(3,*),f(3,*),q(*)
1138
1139      pi = 3.141592654e0
1140      pi43 = 4.0e0*pi/3.0e0
1141      econv = 332.05382e0
1142      fin =  pi43*econv/(xa*yb*zc)
1143
1144      do j=1,3
1145         pole(j) = 0.0e0
1146      end do
1147
1148      do i=1,n
1149         do j=1,3
1150            pole(j) = pole(j) + coo(j,i)*q(i)
1151         end do
1152      end do
1153
1154      poll = 0.0e0
1155      do j=1,3
1156         poll = poll + pole(j)*pole(j)
1157      end do
1158
1159      uin = -0.5e0*fin*poll
1160
1161      cellder(1) = -uin/xa
1162      cellder(3) = -uin/yb
1163      cellder(6) = -uin/zc
1164
1165      do i=1,n
1166         do j=1,3
1167            f(j,i) = f(j,i) - fin*pole(j)*q(i)
1168         end do
1169      end do
1170
1171      return
1172      end
1173
1174      subroutine crionp
1175      parameter (mxion=2000)
1176      implicit real (a-h,p-z),integer (i-n),logical (o)
1177      common /h2oer/   numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat
1178      logical doit
1179
1180      nitmp = 0
1181
1182      do while (nitmp.lt.nion)
1183
1184         doit = .true.
1185         irnd =  numwat*random()
1186
1187         do i=1,nitmp
1188            if (iw(i).eq.irnd) doit = .false.
1189         end do
1190
1191         if (doit) then
1192           nitmp = nitmp + 1
1193           iw(nitmp) = irnd
1194         endif
1195
1196      end do
1197
1198      return
1199      end
1200
1201      real function random()
1202      parameter (mplier=16807,modlus=2147483647,mobymp=127773,
1203     &           momdmp=2836)
1204      common  /seed/jseed,ifrst,nextn
1205c     mseed comes from alloc.c
1206      integer hvlue, lvlue, testv, nextn, mseed
1207
1208      if (ifrst.eq.0) then
1209        jseed = mseed()
1210        nextn = jseed
1211        ifrst = 1
1212      endif
1213
1214      hvlue = nextn / mobymp
1215      lvlue = mod(nextn, mobymp)
1216      testv = mplier*lvlue - momdmp*hvlue
1217
1218      if (testv.gt.0) then
1219        nextn = testv
1220      else
1221        nextn = testv + modlus
1222      endif
1223
1224      random = dble(nextn)/dble(modlus)
1225
1226      return
1227      end
1228
1229