1c
2c $Id$
3c
4
5      function qmmm_geom_distance(r1,r2)
6
7      implicit none
8      double precision r1(3),r2(3)
9      double precision qmmm_geom_distance
10c
11      integer k
12      double precision r
13
14      r=0.0d0
15      do k=1,3
16        r=r+(r1(k)-r2(k))*(r1(k)-r2(k))
17      end do
18      r=sqrt(r)
19
20      qmmm_geom_distance=r
21
22      end
23
24      subroutine qmmm_sort(n,a)
25      implicit none
26      integer n
27      integer a(n)
28c
29c     local variables:
30      integer i
31      integer pass
32      integer sorted
33      integer temp
34      character*32 pname
35
36      pass = 1
37      sorted = 0
38      do while(sorted .eq. 0)
39        sorted = 1
40        do 2 i = 1,n-pass
41          if(a(i) .gt. a(i+1)) then
42            temp = a(i)
43            a(i) = a(i+1)
44            a(i+1) = temp
45            sorted = 0
46          endif
47 2      continue
48        pass = pass +1
49      end do
50      do i=1,n-1
51       if(a(i).eq.a(i+1)) a(i)=-1
52      end do
53
54      return
55
56      end
57
58      function qmmm_map(nr,ar,
59     >                     n,a,map)
60      implicit none
61      integer nr
62      integer ar(nr)
63      integer n
64      integer a(n)
65      integer map(n)
66      logical qmmm_map
67c
68c     local variables:
69      integer i
70      integer ir
71      integer im
72
73      qmmm_map = .true.
74      im = 0
75      do ir=1,nr
76        do i=1,n
77          if(a(i).eq.ar(ir)) then
78            map(i)=ir
79            im = im + 1
80          end if
81        end do
82        if(im.eq.n) return
83      end do
84
85      qmmm_map = .false.
86      return
87      end
88
89      subroutine qmmm_sort_unique(nr,ar,
90     >                     n,a)
91      implicit none
92      integer nr
93      integer ar(nr)
94      integer n
95      integer a(nr)
96c
97c     local variables:
98c
99      integer i
100      integer ir
101      integer im
102
103      n = 0
104      do ir=1,nr
105        do i=1,im
106          if(ar(ir).eq.a(i)) goto 1
107        end do
108        n = n + 1
109        a(n) = ar(ir)
1101       continue
111      end do
112
113      return
114      end
115
116      function qmmm_map1(nr,ar,
117     >                     n,a,nm,map)
118      implicit none
119      integer nr
120      integer ar(nr)
121      integer n
122      integer a(n)
123      integer nm
124      integer map(n)
125      integer qmmm_map1
126c
127c     local variables:
128      integer i
129      integer ir
130
131      nm = 0
132      do ir=1,nr
133        do i=1,n
134          if(a(i).eq.ar(ir)) then
135            map(i)=ir
136            nm = nm + 1
137          end if
138        end do
139        if(nm.eq.n) goto 1
140      end do
141
1421     qmmm_map1 = nm
143
144      return
145      end
146
147      subroutine qmmm_print_bq(nt,c,q)
148      implicit none
149#include "mafdecls.fh"
150#include "errquit.fh"
151#include "qmmm.fh"
152#include "util.fh"
153#include "inp.fh"
154      integer nt
155      double precision c(3,nt)
156      double precision q(nt)
157      character*16 t
158
159      integer i
160
161
162      t = "Bq"
163      do i=1,nt
164       write(6,FMT=9000)
165     >           i,t,c(1,i)*0.529177249d00,
166     >            c(2,i)*0.529177249d00,
167     >            c(3,i)*0.529177249d00, q(i)
168
169      end do
1709000  format(i5,2x,a16,1x,3f15.8,3x,"charge",3x,f15.8)
171
172
173      end
174
175      subroutine qmmm_print_pdb(nt,myname,c,q,t)
176      implicit none
177#include "mafdecls.fh"
178#include "errquit.fh"
179#include "qmmm.fh"
180#include "util.fh"
181#include "inp.fh"
182      integer nt
183      double precision c(3,nt)
184      double precision q(nt)
185      character*16 t(nt)
186      character*(*) myname
187      character*(nw_max_path_len) filename
188
189      integer i
190      integer n
191      integer nf
192      integer ns
193      integer un
194
195      call util_file_prefix(" ",filename)
196      nf=inp_strlen(filename)-1
197      ns=inp_strlen(myname)
198      filename = filename(1:nf)//myname(1:ns)
199c
200      if(.not.qmmm_get_io_unit(un))
201     >   call errquit("cannot get file number",0,0)
202c
203
204      open(unit=un,status="unknown",form="formatted",file=filename)
205
206      do i=1,nt
207       write(un,FMT=9000)
208     >           i,t(i),c(1,i)*0.529177249d00,
209     >            c(2,i)*0.529177249d00,
210     >            c(3,i)*0.529177249d00, q(i)
211
212      end do
2139000  FORMAT("ATOM",T7,I5,T13,A4,T31,F8.3,T39,F8.3,T47,F8.3,T55,F6.2)
214
215      call util_flush(un)
216      close(un)
217
218
219      end
220
221      subroutine qmmm_set_ma_char(nt,myname,t)
222      implicit none
223#include "mafdecls.fh"
224#include "errquit.fh"
225#include "qmmm.fh"
226#include "util.fh"
227#include "inp.fh"
228      integer nt
229      character*16 t(nt)
230      character*(*) myname
231
232      integer i
233
234      do i=1,nt
235       t(i) = myname
236      end do
237
238      end
239
240      subroutine qmmm_print_pdb_bq(nt,myname,c,q)
241      implicit none
242#include "mafdecls.fh"
243#include "errquit.fh"
244#include "qmmm.fh"
245#include "util.fh"
246#include "inp.fh"
247      integer nt
248      double precision c(3,nt)
249      double precision q(nt)
250      character*16 t
251      character*(*) myname
252      character*(nw_max_path_len) filename
253
254      integer i
255      integer n
256      integer nf
257      integer ns
258
259      call util_file_prefix(" ",filename)
260      nf=inp_strlen(filename)-1
261      ns=inp_strlen(myname)
262      filename = filename(1:nf)//myname(1:ns)
263
264      open(unit=46,status="unknown",form="formatted",file=filename)
265
266      t = "Bq"
267      do i=1,nt
268       write(46,FMT=9000)
269     >           i,t,c(1,i)*0.529177249d00,
270     >            c(2,i)*0.529177249d00,
271     >            c(3,i)*0.529177249d00, q(i)
272
273      end do
2749000  FORMAT("ATOM",T7,I5,T13,A4,T31,F8.3,T39,F8.3,T47,F8.3,T55,F6.2)
275
276      call util_flush(46)
277      close(46)
278
279
280      end
281
282      subroutine qmmm_print_pdb_forces(nt,myname,c,f)
283      implicit none
284#include "mafdecls.fh"
285#include "errquit.fh"
286#include "qmmm.fh"
287#include "util.fh"
288#include "inp.fh"
289      integer nt
290      double precision c(3,nt)
291      double precision f(3,nt)
292      character*16 t
293      character*(*) myname
294      character*(nw_max_path_len) filename
295
296      integer i
297      integer n
298      integer nf
299      integer ns
300
301      call util_file_prefix(" ",filename)
302      nf=inp_strlen(filename)-1
303      ns=inp_strlen(myname)
304      filename = filename(1:nf)//myname(1:ns)
305
306      open(unit=46,status="unknown",form="formatted",file=filename)
307
308      t = "Bq"
309      do i=1,nt
310       write(46,FMT=9000)
311     >           i,t,c(1,i)*0.529177249d00,
312     >            c(2,i)*0.529177249d00,
313     >            c(3,i)*0.529177249d00,
314     >            f(1,i),f(2,i),f(3,i)
315
316      end do
3179000  FORMAT("ATOM",T7,I5,T13,A4,T31,F8.3,T39,F8.3,T47,F8.3,T55,3F12.6)
318
319      call util_flush(46)
320      close(46)
321
322
323      end
324
325      subroutine qmmm_print_pdbi(un,nt,ai,c,q,t)
326      implicit none
327#include "mafdecls.fh"
328#include "errquit.fh"
329#include "qmmm.fh"
330#include "util.fh"
331#include "inp.fh"
332      integer un
333      integer nt
334      integer ai(nt)
335      double precision c(3,nt)
336      double precision q(nt)
337      character*16 t(nt)
338
339      integer i
340      integer n
341      integer nf
342      integer ns
343
344      do i=1,nt
345       write(un,FMT=9000)
346     >           ai(i),t(i),c(1,i)*0.529177249d00,
347     >            c(2,i)*0.529177249d00,
348     >            c(3,i)*0.529177249d00, q(i)
349
350      end do
3519000  FORMAT("ATOM",T7,I5,T13,A4,T31,F8.3,T39,F8.3,T47,F8.3,T55,F6.2)
352
353      end
354
355      subroutine qmmm_print_xyz(un,nt,c,t)
356      implicit none
357#include "mafdecls.fh"
358#include "errquit.fh"
359#include "qmmm.fh"
360#include "util.fh"
361#include "inp.fh"
362      integer un
363      integer nt
364      double precision c(3,nt)
365      character*16 t(nt)
366
367      integer i
368      integer n
369      integer nf
370      integer ns
371
372      write(un,FMT=1) nt
373
374      do i=1,nt
375       write(un,FMT=2)
376     >            t(i),c(1,i)*0.529177249d00,
377     >            c(2,i)*0.529177249d00,
378     >            c(3,i)*0.529177249d00
379
380      end do
3811     FORMAT(1X,I9,/, "   ")
3822     FORMAT(1X,A5,6X,3(F12.6,4X))
383
384      end
385
386      subroutine qmmm_print_xyzi(un,nt,ai,c,t)
387      implicit none
388#include "mafdecls.fh"
389#include "errquit.fh"
390#include "qmmm.fh"
391#include "util.fh"
392#include "inp.fh"
393      integer un
394      integer nt
395      integer ai(nt)
396      double precision c(3,nt)
397      character*16 t(nt)
398      character*4 t1
399
400      integer i
401      integer n
402      integer nf
403      integer ns
404
405      write(un,FMT=1) nt
406
407      do i=1,nt
408        t1 = t(i)
409c       write(un,FMT=2)
410       write(un,2)
411     >            t1,c(1,i)*0.529177249d00,
412     >            c(2,i)*0.529177249d00,
413     >            c(3,i)*0.529177249d00,ai(i)
414
415      end do
4161     FORMAT(1X,I9,/, "   ")
4172     FORMAT(1X,A5,6X,3(F12.6,4X),I6)
418
419      end
420
421      function qmmm_get_io_unit(fn)
422
423      implicit none
424      integer fn
425      logical qmmm_get_io_unit
426c
427      integer k
428      logical ostatus
429c
430      do k=80,90
431        INQUIRE(UNIT=k,OPENED=ostatus)
432        ostatus = .not.ostatus
433        if(ostatus)
434     >    INQUIRE(UNIT=k,EXIST=ostatus)
435        if(ostatus) then
436          fn = k
437          qmmm_get_io_unit = .true.
438          return
439        end if
440      end do
441      qmmm_get_io_unit = .false.
442      return
443      end
444
445      function qmmm_file_exist(filename,perm)
446
447      implicit none
448#include "msgids.fh"
449#include "mafdecls.fh"
450#include "global.fh"
451      character*(*) filename
452      logical perm
453      logical qmmm_file_exist
454c
455      character*255 filename0
456      logical ofile
457c
458      filename0 = filename
459      if(perm)
460     >    call util_file_name_resolve(filename0,.false.)
461      if(ga_nodeid().eq.0) then
462        inquire(file=filename0,exist=ofile)
463      end if
464      call ga_brdcst(msg_qmmm_misc, ofile,
465     $     MA_sizeof(MT_INT,1,MT_BYTE),0)
466      call ga_sync()
467c
468      qmmm_file_exist = ofile
469      return
470      end
471
472      subroutine qmmm_interp_xyzi_file(aname1,aname2,aname3,lambda)
473      implicit none
474#include "mafdecls.fh"
475#include "errquit.fh"
476#include "msgids.fh"
477#include "qmmm.fh"
478#include "qmmm_params.fh"
479#include "global.fh"
480#include "inp.fh"
481#include "stdio.fh"
482      character*(*) aname1
483      character*(*) aname2
484      character*(*) aname3
485      double precision lambda
486c     local variables
487      integer ns
488      integer i
489      integer k
490      logical title
491
492
493      integer i_itmp,h_itmp
494      integer i_ctmp,h_ctmp
495      character*32 pname
496      character*30 buf
497      character*72 message
498      character*255 filename1
499      character*255 filename2
500      character*4 tag(1000)
501      double precision r1(3),r2(3)
502      integer ai1,ai2
503
504      integer fn1,fn2
505      integer ns1,ns2
506
507      pname = "qmmm_interp_xyzi_file"
508
509      filename1 = aname1
510      call util_file_name_resolve(filename1,.false.)
511      if(.not.qmmm_get_io_unit(fn1))
512     >       call errquit("cannot get file number",0,0)
513c
514      message = "opening filename1 for reading"
515      open(fn1,file=filename1,form='formatted',status='old',
516     $          err=133)
517c
518      filename2 = aname2
519      call util_file_name_resolve(filename2,.false.)
520      if(.not.qmmm_get_io_unit(fn2))
521     >       call errquit("cannot get file number",0,0)
522c
523      message = "opening filename2 for reading"
524      open(fn2,file=filename2,form='formatted',status='old',
525     $          err=133)
526
527c     get number of atoms
528      read(fn1,*) ns1
529      read(fn2,*) ns2
530      if(ns1.ne.ns2)
531     >  call errquit( pname//'different number of atoms',0,0)
532      ns = ns1
533      if(ns.gt.1000)
534     > call errquit(pname//"increase memory for atom tag",ns,0)
535
536      if(.not.ma_push_get(mt_int,ns,'itmp',h_itmp,i_itmp))
537     + call errquit( pname//'Failed to allocate memory for itmp1',
538     + ns, MA_ERR)
539      call ifill(ns,0,int_mb(i_itmp),1)
540
541      if(.not.ma_push_get(mt_dbl,3*ns,'ctmp',h_ctmp,i_ctmp))
542     + call errquit( pname//'Failed to allocate memory for ctmp1',
543     + 3*ns, MA_ERR)
544      call dfill(3*ns,0.0d0,dbl_mb(i_ctmp),1)
545
546c
547c     forward to cooordinates
548c     -----------------------
549      rewind(fn1)
550      title =.false.
551      do i=1,ns+2
552        read(fn1,*,end=1) buf
553      end do
554      title = .true.
5551     continue
556      rewind(fn1)
557      read(fn1,*) buf
558      if(title)
559     +   read(fn1,*) buf
560
561      rewind(fn2)
562      title =.false.
563      do i=1,ns+2
564        read(fn2,*,end=2) buf
565      end do
566      title = .true.
5672     continue
568      rewind(fn2)
569      read(fn2,*) buf
570      if(title)
571     +   read(fn2,*) buf
572
573      do i=1,ns
574       read(fn1,*) tag(i), (r1(k),k=1,3),
575     +             ai1
576       read(fn2,*) buf, (r2(k),k=1,3),
577     +             ai2
578
579       if(ai2.ne.ai1)
580     >   call errquit( pname//'different global index',0,0)
581       int_mb(i_itmp+i-1) = ai1
582       do k=1,3
583        dbl_mb(i_ctmp+3*(i-1)+k-1) = r1(k)+lambda*(r2(k)-r1(k))
584       end do
585      end do
586c
587c
588      filename1 = aname3
589      call util_file_name_resolve(filename1,.false.)
590c
591      message = "opening filename1 for writing"
592      open(fn1,file=filename1,form='formatted',status='unknown',
593     $          err=133)
594
595      write(fn1,3) ns
596      do i=1,ns
597       write(fn1,4) tag(i), (dbl_mb(i_ctmp+3*(i-1)+k-1),k=1,3),
598     +             int_mb(i_itmp+i-1)
599
600      end do
6013     FORMAT(1X,I9,/, "   ")
6024     FORMAT(1X,A5,6X,3(F12.6,4X),I6)
603      close(fn1)
604      close(fn2)
605
606c
607      if(.not.ma_pop_stack(h_ctmp))
608     & call errquit('qmmm:
609     >              Failed to deallocate stack c_tmp',ns,
610     &       MA_ERR)
611
612      if(.not.ma_pop_stack(h_itmp))
613     & call errquit('qmmm:
614     >              Failed to deallocate stack i_itmp',ns,
615     &       MA_ERR)
616
617
618      return
619
620 133  call errquit(pname//'error '//message,0, 0)
621
622      end
623
624      subroutine qmmm_interp_esp_file(aname1,aname2,aname3,lambda)
625      implicit none
626#include "mafdecls.fh"
627#include "errquit.fh"
628#include "msgids.fh"
629#include "qmmm.fh"
630#include "qmmm_params.fh"
631#include "global.fh"
632#include "inp.fh"
633#include "stdio.fh"
634      character*(*) aname1
635      character*(*) aname2
636      character*(*) aname3
637      double precision lambda
638c     local variables
639      integer ns
640      integer i
641      integer k
642      logical title
643
644
645      integer i_itmp,h_itmp
646      integer i_ctmp,h_ctmp
647      character*32 pname
648      character*30 buf
649      character*72 message
650      character*255 filename1
651      character*255 filename2
652      double precision q1,q2,q3(1000)
653      integer ai1,ai2,ai3(1000)
654
655      integer fn1,fn2
656      integer ns1,ns2
657
658      pname = "qmmm_interp_esp_file"
659
660      filename1 = aname1
661      call util_file_name_resolve(filename1,.false.)
662      if(.not.qmmm_get_io_unit(fn1))
663     >       call errquit("cannot get file number",0,0)
664c
665      message = "opening filename1 for reading"
666      open(fn1,file=filename1,form='formatted',status='old',
667     $          err=133)
668c
669      filename2 = aname2
670      call util_file_name_resolve(filename2,.false.)
671      if(.not.qmmm_get_io_unit(fn2))
672     >       call errquit("cannot get file number",0,0)
673c
674      message = "opening filename2 for reading"
675      open(fn2,file=filename2,form='formatted',status='old',
676     $          err=133)
677
678
679      ns = 0
680      do i=1,1000
681       read(fn1,*,end=1) ai1,q1
682       read(fn2,*,end=1) ai2,q2
683
684       ns = ns + 1
685       if(ai2.ne.ai1)
686     >   call errquit( pname//'different global index',0,0)
687       if(ns.gt.1000)
688     > call errquit(pname//"increase memory for atom tag",ns,0)
689       ai3(i) = ai1
690       q3(i) = q1+lambda*(q2-q1)
691      end do
692c
693 1    continue
694c
695c
696      filename1 = aname3
697      call util_file_name_resolve(filename1,.false.)
698c
699      message = "opening filename1 for writing"
700      open(fn1,file=filename1,form='formatted',status='unknown',
701     $          err=133)
702
703      do i=1,ns
704       write(fn1,*) ai3(i),q3(i)
705      end do
706
707      close(fn1)
708      close(fn2)
709
710
711      return
712
713 133  call errquit(pname//'error '//message,0, 0)
714
715      end
716