1      subroutine qmmm_data_init(irtdb)
2      implicit none
3c
4#include "mafdecls.fh"
5#include "errquit.fh"
6#include "qmmm_data.fh"
7#include "qmmm_params.fh"
8#include "rtdb.fh"
9#include "global.fh"
10#include "inp.fh"
11
12      integer irtdb
13c
14      character*32 pname
15      logical qmmm_h_link
16      external qmmm_h_link
17
18
19      pname = "qmmm_data_init: "
20
21      call mm_get_tot_nqm(nqm)
22      if(ga_nodeid().eq.0)
23     >  write(*,*) "number of quantum atoms",nqm
24c
25c     allocate data arrays
26c     --------------------
27      call qmmm_data_alloc()
28c
29c     establish persistent mapping between mm and qm atoms
30c     first go pure quantum atoms and then links
31c     ----------------------------------------------------
32      call mm_get_solute_quant_ind(nqm,int_mb(i_ai))
33      if(.not.rtdb_get(irtdb,'charge',mt_dbl,1,qcharge)) then
34        qcharge = 0.0d0
35      end if
36
37      end
38
39      subroutine qmmm_data_rdinp(irtdb)
40      implicit none
41c
42#include "mafdecls.fh"
43#include "errquit.fh"
44#include "qmmm_data.fh"
45#include "qmmm_link_data.fh"
46#include "qmmm_params.fh"
47#include "rtdb.fh"
48#include "global.fh"
49#include "inp.fh"
50#include "util.fh"
51
52      integer irtdb
53c
54      character*84 tag
55      character*32 pname
56      character*32 optimization
57      character*32 link_ecp
58      double precision  eatoms
59      integer linkatm,bq_exclude
60      character*30 operation
61      character*30 region(3)
62      integer nregion
63      logical ignore
64      integer i
65
66      pname = "qmmm_data_init: "
67
68      master = ga_nodeid().eq.0
69c
70      oprint_debug = util_print("debug_info", print_debug)
71      oprint_debug = oprint_debug .and. (ga_nodeid().eq.0)
72c
73      oprint_low = util_print("low", print_low)
74c
75      eatoms = 0.0d0
76      linkatm=qmmm_hydrogen
77      bq_exclude=qmmm_no_atoms
78      link_ecp = "auto"
79c
80c     initialize rtdb parameters
81      if(master)
82     > call util_print_centered(6,
83     >    "QM/MM Interface Parameters",32,.true.)
84c
85c     get the operation
86c     -----------------
87      tag='task:operation'
88      if (.not. rtdb_cget(irtdb, 'task:operation', 1, operation))
89     $     operation = 'energy'
90      if (master)
91     >  write(6,22) "operation",tag,operation
92c
93      tag = "qmmm:uqmatm"
94      if (.not.rtdb_get(irtdb,"qmmm:uqmatm",mt_dbl,1,eatoms)) then
95        if (.not.rtdb_put(irtdb,"qmmm:uqmatm",mt_dbl,1,eatoms))
96     >      call errquit(pname//'failed to store eatoms',0,
97     >       RTDB_ERR)
98      end if
99      if(master) then
100       if(eatoms.eq.0) then
101         write(6,19) "reference energy",tag,eatoms,
102     >     "  <--- Warning zero value is not advisable !!!"
103       else
104         write(6,20) "reference energy",tag,eatoms
105       end if
106      end if
107
108c
109c     default value for bqzone from md
110      if(.not.rtdb_get(irtdb,'md:rqmmm',mt_dbl,1,bqzone))
111     + call errquit(pname//'failed to get bqzone from mm',0)
112c     convert from nm to angstrom
113      bqzone = bqzone*cnm2ang
114c
115      tag = "qmmm:bqzone"
116      if (rtdb_get(irtdb,tag,mt_dbl,1,bqzone)) then
117        if (.not.rtdb_put(irtdb,'md:rqmmm',mt_dbl,1,bqzone/cnm2ang))
118     >      call errquit(pname//'failed to store '//tag,0,
119     >       RTDB_ERR)
120      end if
121      if(master)
122     > write(6,20) "bqzone radius",tag,bqzone
123
124
125      tag ="qmmm:bq_exclude"
126      if (.not.rtdb_get(irtdb,tag,mt_int,1,bq_exclude)) then
127        if (.not.rtdb_put(irtdb,tag,mt_int,1,bq_exclude))
128     >      call errquit(pname//'failed to store bq_exclude',0,
129     >       RTDB_ERR)
130
131       end if
132      if(master) then
133        if(bq_exclude.eq.qmmm_linkbond_H) then
134          write(6,22) "excluded bq's",tag,
135     >                 "hydrogens bonded to link atoms"
136        else if(bq_exclude.eq.qmmm_linkbond) then
137          write(6,22) "excluded bq's",tag,"all atoms bonded to links"
138        else if(bq_exclude.eq.qmmm_all_atoms) then
139          write(6,22) "excluded bq's",tag,"all"
140        else if(bq_exclude.eq.qmmm_no_atoms) then
141          write(6,22) "excluded bq's",tag,"none"
142        else
143          call errquit(pname//'invalid bq_exclude',0,RTDB_ERR)
144        end if
145      end if
146
147      tag ="qmmm:linkatm"
148      if (.not.rtdb_get(irtdb,tag,mt_int,1,linkatm)) then
149        if (.not.rtdb_put(irtdb,tag,mt_int,1,linkatm))
150     >      call errquit(pname//'failed to store'//tag,0,
151     >       RTDB_ERR)
152
153       end if
154       if(master) then
155         if(linkatm.eq.qmmm_hydrogen) then
156           write(6,22) "link atom type",tag,"hydrogens"
157         else if(linkatm.eq.qmmm_halogen) then
158           write(6,22) "link atom type",tag,"halogens"
159         else
160           call errquit(pname//'invalid link atom type',0,RTDB_ERR)
161         end if
162       end if
163
164      optimization = "bfgs"
165      tag ="qmmm:optimization"
166      if (.not.rtdb_cget(irtdb,tag,1,optimization)) then
167        if (.not.rtdb_cput(irtdb,tag,1,optimization))
168     >      call errquit(pname//'failed to store'//tag,0,
169     >       RTDB_ERR)
170
171      end if
172      if (master)
173     >  write(6,22) "optimization method",tag,optimization
174
175c
176c     region definitions
177c     ------------------
178      tag ="qmmm:region"
179      if (.not.rtdb_get(irtdb,tag(1:inp_strlen(tag))//"_n",
180     >                 mt_int,1,nregion)) then
181c        if(operation.ne."energy" .and.
182c     >     operation.ne."gradient" .and.
183c     >     operation.ne."property" .and.
184c     >     operation.ne."fep" .and.
185c     >     operation.ne."neb" )
186c     >     call errquit(pname//'cannot find '//tag,0,
187c     >       RTDB_ERR)
188        nregion = 0
189      end if
190      if(nregion.ne.0) then
191        if (.not.rtdb_cget(irtdb,tag,nregion,region))
192     >      call errquit(pname//"failed to get"//tag,0,RTDB_ERR)
193      end if
194      if (.not.rtdb_put(irtdb,tag(1:inp_strlen(tag))//"_n",
195     >                  mt_int,1,nregion))
196     >      call errquit(pname//'failed to store'//tag,0,
197     >                   RTDB_ERR)
198
199      if (master) then
200        tag = " "
201        do i=1,nregion
202          write(tag,'("qmmm:region ",I1)') i
203          write(6,22) tag,region(i)
204        end do
205      end if
206
207      tag ="qmmm:link_ecp"
208      if (.not.rtdb_cget(irtdb,tag,1,link_ecp)) then
209       if (.not.rtdb_cput(irtdb,tag,1,link_ecp))
210     >      call errquit(pname//'failed to store'//tag,0,
211     >       RTDB_ERR)
212
213       else if (link_ecp .ne. "auto" .and.
214     >          link_ecp .ne. "user" ) then
215            call errquit(pname//'unknown value'//tag,0,
216     >       RTDB_ERR)
217
218       end if
219      if (master)
220     >  write(6,22) "ecp on link atoms",tag,link_ecp
221
222      qmmm_interface='qm'
223      tag ="qmmm:interface"
224      ignore = rtdb_cget(irtdb,tag,1,qmmm_interface)
225      if (.not.rtdb_cget(irtdb,tag,1,qmmm_interface)) then
226         qmmm_interface = "qm"
227      end if
228      if (master)
229     >  write(6,22) "interface api",tag,qmmm_interface
230
231      tag = "qmmm:bq_dynamical"
232      if (.not.rtdb_get(irtdb,tag,mt_log,1,bq_dynamical))
233     >  bq_dynamical = .false.
234
235      if (master)
236     >  write(6,23)
237
238
23919    FORMAT(1X,A,T24,A,T46,F12.6,A)
24020    FORMAT(1X,A,T24,A,T46,F12.6)
24121    FORMAT(1X,A,T24,A,T46,L3)
24222    FORMAT(1X,A,T24,A,T46,A)
24323    FORMAT(1X,54("-"),//)
244
245       if (.not.rtdb_get(irtdb,'qmmm:linkatm',mt_int,1,link_atom))
246     + call errquit('qmmm_data_init: qmmm:linkatm',link_atom,
247     &       RTDB_ERR)
248
249      nqm   = -1
250
251      end
252
253      subroutine qmmm_data_alloc()
254      implicit none
255c
256#include "mafdecls.fh"
257#include "errquit.fh"
258#include "qmmm_data.fh"
259c
260c
261c     indexing array
262c     -------------
263      if(.not.ma_alloc_get(MT_INT, nqm, 'qmmm index array',
264     &      h_ai, i_ai) ) call errquit(
265     &      'qmmm_data_alloc: unable to allocate heap space',
266     &      nqm, MA_ERR)
267      call ifill(nqm,-1,int_mb(i_ai),1)
268
269      end
270
271      subroutine qmmm_data_release()
272      implicit none
273#include "mafdecls.fh"
274#include "errquit.fh"
275#include "qmmm_data.fh"
276
277      if(.not.ma_free_heap(h_ai))
278     & call errquit('qmmm ai: Failed to deallocate heap',nqm,
279     &       MA_ERR)
280
281
282       return
283      end
284
285      function qmmm_get_nqm()
286      implicit none
287#include "qmmm_data.fh"
288
289      integer qmmm_get_nqm
290
291      qmmm_get_nqm = nqm
292
293      end
294
295      function qmmm_get_i_ai()
296      implicit none
297#include "qmmm_data.fh"
298      integer qmmm_get_i_ai
299
300      qmmm_get_i_ai = i_ai
301
302      end
303
304      function qmmm_master()
305      implicit none
306#include "qmmm_data.fh"
307#include "qmmm_params.fh"
308
309      logical qmmm_master
310
311      qmmm_master = master
312
313      end
314
315      function qmmm_bq_dynamical()
316      implicit none
317#include "qmmm_data.fh"
318#include "qmmm_params.fh"
319
320      logical qmmm_bq_dynamical
321
322      qmmm_bq_dynamical = bq_dynamical
323
324      end
325
326      function qmmm_print_default(name)
327      implicit none
328#include "qmmm_data.fh"
329#include "util.fh"
330#include "global.fh"
331
332      logical qmmm_print_default
333      character*(*) name
334c
335      logical oprint
336
337      oprint = util_print(name, print_default)
338
339      qmmm_print_default = oprint .and. (ga_nodeid().eq.0)
340
341      end
342
343      function qmmm_print_debug()
344      implicit none
345#include "qmmm_data.fh"
346#include "util.fh"
347#include "global.fh"
348
349      logical qmmm_print_debug
350c
351
352c      oprint = util_print("debug_info", print_debug)
353c      qmmm_print_debug = oprint .and. (ga_nodeid().eq.0)
354
355       qmmm_print_debug = oprint_debug
356      end
357
358      function qmmm_print_low()
359      implicit none
360#include "qmmm_data.fh"
361#include "util.fh"
362#include "global.fh"
363
364      logical qmmm_print_low
365c
366
367       qmmm_print_low = oprint_low
368      end
369
370      function qmmm_get_bqzone()
371      implicit none
372#include "qmmm_data.fh"
373      double precision qmmm_get_bqzone
374
375      qmmm_get_bqzone = bqzone
376
377      end
378
379      function qmmm_get_interface()
380      implicit none
381#include "qmmm_data.fh"
382      character*32 qmmm_get_interface
383
384      qmmm_get_interface = qmmm_interface
385
386      end
387
388      subroutine qmmm_set_interface(a_interface)
389      implicit none
390#include "qmmm_data.fh"
391      character*(*) a_interface
392
393      qmmm_interface = a_interface
394
395      end
396
397c $Id$
398