1      subroutine mm_vdw_coords_init()
2      implicit none
3#include "util.fh"
4#include "errquit.fh"
5#include "inp.fh"
6#include "stdio.fh"
7#include "rtdb.fh"
8#include "mafdecls.fh"
9#include "mm_coords_data.fh"
10#include "mm_vdw_data.fh"
11#include "mm_vdw_coords_data.fh"
12      integer i,j
13      integer nvtot0,nvqm0,nvmm0
14      integer nv14tot0,nv14qm0,nv14mm0
15      character*30 pname
16
17      pname="mm_vdw_coords_init"
18c     write(*,*) pname
19
20      if(.not.ma_alloc_get(mt_log,nqm,"add qm to vdw coords list",
21     &   h_advqm, i_advqm)) goto 911
22
23      if(.not.ma_alloc_get(mt_log,nmm,"add mm to vdw coords list",
24     &   h_advmm, i_advmm)) goto 911
25
26      if(.not.ma_alloc_get(mt_log,nqm,"add qm to vdw14 coords list",
27     &   h_advqm14, i_advqm14)) goto 911
28
29      if(.not.ma_alloc_get(mt_log,nmm,"add mm to vdw14 coords list",
30     &   h_advmm14, i_advmm14)) goto 911
31
32      call mm_vdw_coords_map(nvqm0,nqm,log_mb(i_advqm),
33     &   int_mb(i_iqm),nvdw,int_mb(i_ivdw),int_mb(i_jvdw))
34
35      call mm_vdw_coords_map(nvmm0,nmm,log_mb(i_advmm),
36     &   int_mb(i_imm),nvdw,int_mb(i_ivdw),int_mb(i_jvdw))
37
38      call mm_vdw_coords_map(nv14qm0,nqm,log_mb(i_advqm14),
39     &   int_mb(i_iqm),nvdw14,int_mb(i_ivdw14),int_mb(i_jvdw14))
40
41      call mm_vdw_coords_map(nv14mm0,nmm,log_mb(i_advmm14),
42     &   int_mb(i_imm),nvdw14,int_mb(i_ivdw14),int_mb(i_jvdw14))
43
44      call mm_vdw_coords_allocate(nvqm0,nvmm0)
45      call mm_vdw14_coords_allocate(nv14qm0,nv14mm0)
46
47      call mm_vdw_coords_load(nvtot,nvqm,nvmm,dbl_mb(i_rvdw),
48     &     int_mb(i_icvdw),log_mb(i_lvqm),
49     &     log_mb(i_advqm),log_mb(i_advmm),.false.)
50
51      call mm_vdw_coords_load(nv14tot,nv14qm,nv14mm,dbl_mb(i_rvdw14),
52     &     int_mb(i_icvdw14),log_mb(i_lv14qm),
53     &     log_mb(i_advqm14),log_mb(i_advmm14),.true.)
54
55      if(.not.ma_free_heap(h_advqm))
56     & call errquit('mm:
57     &              Failed to deallocate stack advqm',nqm,
58     &       MA_ERR)
59
60      if(.not.ma_free_heap(h_advmm))
61     & call errquit('mm:
62     &              Failed to deallocate stack advmm',nmm,
63     &       MA_ERR)
64
65      if(.not.ma_free_heap(h_advqm14))
66     & call errquit('mm:
67     &              Failed to deallocate stack advqm14',nqm,
68     &       MA_ERR)
69
70      if(.not.ma_free_heap(h_advmm14))
71     & call errquit('mm:
72     &              Failed to deallocate stack advqm14',nmm,
73     &       MA_ERR)
74
75      return
76911   call errquit("error "//trim(pname),0,
77     >        -1)
78
79      end
80
81      subroutine mm_vdw_coords_map(nvt,nt,addtyp,
82     &           indx_typ,nij,i_i,i_j)
83      implicit none
84#include "util.fh"
85#include "errquit.fh"
86#include "inp.fh"
87#include "stdio.fh"
88#include "rtdb.fh"
89#include "mafdecls.fh"
90      integer nvt
91      integer nt
92      logical addtyp(nt)
93      integer indx_typ(nt)
94      integer nij
95      integer i_i(nij),i_j(nij)
96
97      integer i,j
98      integer indx_i, indx_j
99
100      nvt = 0
101      addtyp(1:nt) = .false.
102
103      do i=1,nt
104        do j=1,nij
105          indx_i = i_i(j)
106          indx_j = i_j(j)
107          if(indx_i.eq.indx_typ(i)) then
108            if(.not.addtyp(i)) then
109              addtyp(i) = .true.
110              nvt = nvt + 1
111            end if
112          end if
113
114          if(indx_j.eq.indx_typ(i)) then
115            if(.not.addtyp(i)) then
116              addtyp(i) = .true.
117              nvt = nvt + 1
118            end if
119          end if
120        end do
121      end do
122
123      end
124
125      subroutine mm_vdw_coords_allocate(nvqm0,nvmm0)
126      implicit none
127#include "errquit.fh"
128#include "stdio.fh"
129#include "mafdecls.fh"
130#include "mm_vdw_coords_data.fh"
131      integer nvqm0, nvmm0
132
133      character*180 message
134      character*30 pname
135
136      integer nvtot0
137
138      pname = "mm_vdw_coords_allocate"
139      nvtot0 = nvqm0 + nvmm0
140      if(nvtot0.ne.nvtot) then
141        call mm_vdw_coords_end()
142        if(.not.ma_alloc_get(mt_dbl,3*nvtot0,
143     &                        "mm vdw coords",
144     &                        h_rvdw,i_rvdw)) goto 911
145        if(.not.ma_alloc_get(mt_int,nvtot0,
146     &                        "mm vdw indices",
147     &                        h_icvdw,i_icvdw)) goto 911
148        if(.not.ma_alloc_get(mt_log,nvtot0,
149     &                        "mm vdw quantum flags",
150     &                        h_lvqm,i_lvqm)) goto 911
151      end if
152      nvtot = nvtot0
153      nvqm = nvqm0
154      nvmm = nvmm0
155      call dfill(3*nvtot,0.0d0,dbl_mb(i_rvdw),1)
156      call ifill(nvtot,0.0,int_mb(i_icvdw),1)
157
158      return
159911   call errquit("error "//trim(pname),0,-1)
160      end
161
162      subroutine mm_vdw14_coords_allocate(nv14qm0,nv14mm0)
163      implicit none
164#include "errquit.fh"
165#include "stdio.fh"
166#include "mafdecls.fh"
167#include "mm_vdw_coords_data.fh"
168      integer nv14qm0, nv14mm0
169
170      character*180 message
171      character*30 pname
172
173      integer nv14tot0
174
175      pname = "mm_vdw14_coords_allocate"
176      nv14tot0 = nv14qm0 + nv14mm0
177      if(nv14tot0.ne.nv14tot) then
178        call mm_vdw14_coords_end()
179        if(.not.ma_alloc_get(mt_dbl,3*nv14tot0,
180     &                        "mm vdw14 coords",
181     &                        h_rvdw14,i_rvdw14)) goto 911
182        if(.not.ma_alloc_get(mt_int,nv14tot0,
183     &                        "mm vdw14 indices",
184     &                        h_icvdw14,i_icvdw14)) goto 911
185        if(.not.ma_alloc_get(mt_dbl,nv14tot0,
186     &                        "mm vdw14 charges",
187     &                        h_chgmm14,i_chgmm14)) goto 911
188        if(.not.ma_alloc_get(mt_log,nv14tot0,
189     &                        "mm vdw14 quantum flags",
190     &                        h_lv14qm,i_lv14qm)) goto 911
191      end if
192      nv14tot = nv14tot0
193      nv14qm = nv14qm0
194      nv14mm = nv14mm0
195      call dfill(3*nv14tot,0.0d0,dbl_mb(i_rvdw14),1)
196      call dfill(nv14tot,0.0d0,dbl_mb(i_chgmm14),1)
197      call ifill(nv14tot,0.0,int_mb(i_icvdw14),1)
198
199      return
200911   call errquit("error "//trim(pname),0,-1)
201      end
202
203      subroutine mm_vdw_coords_end()
204      implicit none
205#include "util.fh"
206#include "errquit.fh"
207#include "inp.fh"
208#include "stdio.fh"
209#include "rtdb.fh"
210#include "mafdecls.fh"
211#include "mm_vdw_coords_data.fh"
212
213      character*30 pname
214      pname = "mm_vdw_coords_end"
215
216      if(nvtot.gt.0) then
217          if (.not.ma_free_heap(h_rvdw))   goto 911
218          if (.not.ma_free_heap(h_icvdw))   goto 911
219          if (.not.ma_free_heap(h_lvqm))   goto 911
220          nvtot = 0
221          nvqm = 0
222          nvmm = 0
223      end if
224
225      return
226911   call errquit("error "//trim(pname),0,-1)
227
228      end
229
230      subroutine mm_vdw14_coords_end()
231      implicit none
232#include "util.fh"
233#include "errquit.fh"
234#include "inp.fh"
235#include "stdio.fh"
236#include "rtdb.fh"
237#include "mafdecls.fh"
238#include "mm_vdw_coords_data.fh"
239
240      character*30 pname
241      pname = "mm_vdw14_coords_end"
242
243      if(nv14tot.gt.0) then
244          if (.not.ma_free_heap(h_rvdw14))   goto 911
245          if (.not.ma_free_heap(h_icvdw14))   goto 911
246          if (.not.ma_free_heap(h_chgmm14))   goto 911
247          if (.not.ma_free_heap(h_lv14qm))   goto 911
248          nv14tot = 0
249          nv14qm = 0
250          nv14mm = 0
251      end if
252
253      return
254911   call errquit("error "//trim(pname),0,-1)
255
256      end
257
258      subroutine mm_vdw_coords_deallocate()
259      implicit none
260#include "errquit.fh"
261#include "stdio.fh"
262#include "mafdecls.fh"
263        call mm_vdw_coords_end()
264        call mm_vdw14_coords_end()
265
266      end
267
268      subroutine mm_vdw_coords_load(nvt,nvq,nvm,c,ic,lvqm,
269     +                              advqm,advmm,load_mmchg)
270      implicit none
271#include "util.fh"
272#include "errquit.fh"
273#include "inp.fh"
274#include "stdio.fh"
275#include "rtdb.fh"
276#include "mafdecls.fh"
277#include "mm_coords_data.fh"
278#include "mm_vdw_coords_data.fh"
279      integer nvt,nvq,nvm
280      double precision c(3*nvt)
281      integer ic(nvt)
282      logical lvqm(nvt)
283      logical advqm(nqm)
284      logical advmm(nmm)
285      logical load_mmchg
286
287      character*30 pname
288      integer n
289      integer i
290      integer indx_qm, indx_mm
291      logical addatom
292      double precision coord(3)
293
294      pname = "mm_vdw_coords_load"
295
296      lvqm(1:nvt) = .false.
297      n = 0
298
299      do i=1,nqm
300        addatom = advqm(i)
301        indx_qm = int_mb(i_iqm+i-1)
302        coord(1) = dbl_mb(i_rqm+3*(i-1))
303        coord(2) = dbl_mb(i_rqm+3*(i-1)+1)
304        coord(3) = dbl_mb(i_rqm+3*(i-1)+2)
305        if(addatom) then
306          n = n + 1
307          lvqm(n) = .true.
308          ic(n) = indx_qm
309          c(3*(n-1)+1)   = coord(1)
310          c(3*(n-1)+2)   = coord(2)
311          c(3*(n-1)+3)   = coord(3)
312        end if
313        if(n == nvq) exit
314      end do
315
316      n = 0
317
318      do i=1,nmm
319        addatom = advmm(i)
320        indx_mm = int_mb(i_imm+i-1)
321        coord(1) = dbl_mb(i_rmm+3*(i-1))
322        coord(2) = dbl_mb(i_rmm+3*(i-1)+1)
323        coord(3) = dbl_mb(i_rmm+3*(i-1)+2)
324        if(addatom) then
325          n = n + 1
326          ic(nvq+n) = indx_mm
327          c(3*(nvq+n-1)+1)   = coord(1)
328          c(3*(nvq+n-1)+2)   = coord(2)
329          c(3*(nvq+n-1)+3)   = coord(3)
330          if(load_mmchg) then
331            dbl_mb(i_chgmm14+nvq+n-1) = dbl_mb(i_chgmm+i-1)
332          end if
333        end if
334        if(n == nvm) exit
335      end do
336
337c     call mm_vdw_coords_test(nvt,c,ic)
338      return
339911   call errquit("error "//trim(pname),0,
340     >        -1)
341      return
342      end
343
344      subroutine mm_vdw_qmcoords_load(rtdb)
345      implicit none
346#include "util.fh"
347#include "errquit.fh"
348#include "inp.fh"
349#include "stdio.fh"
350#include "rtdb.fh"
351#include "mafdecls.fh"
352#include "geom.fh"
353#include "mm_geom_data.fh"
354#include "mm_coords_data.fh"
355#include "mm_vdw_coords_data.fh"
356      integer rtdb
357
358      integer n,fn
359      integer geom
360      integer nat
361      logical status
362      character*30 message
363      character*30 pname
364      double precision scale
365
366c     local variables
367      character*16 t(nqm)
368      integer i, j
369      integer iqm, ivqm
370      integer i_ctmp, h_ctmp
371
372      pname = "mm_vdw_qmcoords_load"
373c     write(*,*) pname
374
375c     load geometry
376c     -------------
377      if (.not. geom_create(geom, 'geometry'))
378     &     call errquit('cons_load_geom',0, GEOM_ERR)
379      if (.not. geom_rtdb_load(rtdb, geom, 'geometry'))
380     &     call errquit('cons_load_geom',0, RTDB_ERR)
381
382c     get cart coordinates
383c     --------------------
384      status=geom_ncent(geom,nat)
385
386      if(.not.ma_push_get(mt_dbl,3*nat,'ctmp',h_ctmp,i_ctmp))
387     & call errquit( pname//'Failed to allocate memory for ctmp',
388     & 3*nat, MA_ERR)
389
390      if(.not. geom_cart_coords_get(geom,dbl_mb(i_ctmp)))
391     &    call errquit("mm:geom_cart_coords_get",0,0)
392
393      call util_convert_units("au","angstrom",scale)
394      call dscal(3*nact, scale,dbl_mb(i_ctmp),1)
395
396      do i=1,nact
397        iqm = int_mb(i_iact+i-1)
398        do j=1,nvtot
399          ivqm = int_mb(i_icvdw+j-1)
400          if(ivqm.eq.iqm) then
401            dbl_mb(i_rvdw+3*(j-1)) = dbl_mb(i_ctmp+3*(i-1))
402            dbl_mb(i_rvdw+3*(j-1)+1) = dbl_mb(i_ctmp+3*(i-1)+1)
403            dbl_mb(i_rvdw+3*(j-1)+2) = dbl_mb(i_ctmp+3*(i-1)+2)
404          end if
405        end do
406      end do
407
408      if(.not.ma_pop_stack(h_ctmp))
409     & call errquit('mm:
410     &              Failed to deallocate stack c_tmp',nat,
411     &       MA_ERR)
412
413      if(.not.geom_destroy(geom))
414     &    goto 911
415
416      return
417
418911   call errquit("error "//trim(pname),0,-1)
419
420      end
421
422      subroutine mm_vdw14_qmcoords_load(rtdb)
423      implicit none
424#include "util.fh"
425#include "errquit.fh"
426#include "inp.fh"
427#include "stdio.fh"
428#include "rtdb.fh"
429#include "mafdecls.fh"
430#include "geom.fh"
431#include "mm_geom_data.fh"
432#include "mm_coords_data.fh"
433#include "mm_vdw_coords_data.fh"
434      integer rtdb
435
436      integer n,fn
437      integer geom
438      integer nat
439      logical status
440      character*30 message
441      character*30 pname
442      double precision scale
443
444c     local variables
445      character*16 t(nqm)
446      integer i, j
447      integer iqm, ivqm
448      integer i_ctmp, h_ctmp
449
450      pname = "mm_vdw14_qmcoords_load"
451c     write(*,*) pname
452
453c     load geometry
454c     -------------
455      if (.not. geom_create(geom, 'geometry'))
456     &     call errquit('cons_load_geom',0, GEOM_ERR)
457      if (.not. geom_rtdb_load(rtdb, geom, 'geometry'))
458     &     call errquit('cons_load_geom',0, RTDB_ERR)
459
460c     get cart coordinates
461c     --------------------
462      status=geom_ncent(geom,nat)
463
464      if(.not.ma_push_get(mt_dbl,3*nat,'ctmp',h_ctmp,i_ctmp))
465     & call errquit( pname//'Failed to allocate memory for ctmp',
466     & 3*nat, MA_ERR)
467
468      if(.not. geom_cart_coords_get(geom,dbl_mb(i_ctmp)))
469     &    call errquit("mm:geom_cart_coords_get",0,0)
470
471      call util_convert_units("au","angstrom",scale)
472      call dscal(3*nact, scale,dbl_mb(i_ctmp),1)
473
474      do i=1,nact
475        iqm = int_mb(i_iact+i-1)
476        do j=1,nv14tot
477          ivqm = int_mb(i_icvdw14+j-1)
478          if(ivqm.eq.iqm) then
479            dbl_mb(i_rvdw14+3*(j-1)) = dbl_mb(i_ctmp+3*(i-1))
480            dbl_mb(i_rvdw14+3*(j-1)+1) = dbl_mb(i_ctmp+3*(i-1)+1)
481            dbl_mb(i_rvdw14+3*(j-1)+2) = dbl_mb(i_ctmp+3*(i-1)+2)
482          end if
483        end do
484      end do
485
486      if(.not.ma_pop_stack(h_ctmp))
487     & call errquit('mm:
488     &              Failed to deallocate stack c_tmp',nat,
489     &       MA_ERR)
490
491      if(.not.geom_destroy(geom))
492     &    goto 911
493
494      return
495
496911   call errquit("error "//trim(pname),0,-1)
497
498      end
499
500
501      subroutine mm_vdw_coords_test(nvt,c,ic)
502      implicit none
503#include "mm_vdw_coords_data.fh"
504#include "mafdecls.fh"
505      integer nvt
506      double precision c(nvt)
507      integer ic(nvt)
508
509      integer i,j
510      do i=1,nvt
511         write(6,'(1X,3(1X,F10.6),I5)')
512     $        (c(3*(i-1)+j),j=1,3),
513     $        ic(i)
514      end do
515
516      end
517