1      subroutine qmmm_check_forces1(rtdb)
2      implicit none
3#include "errquit.fh"
4#include "global.fh"
5#include "mafdecls.fh"
6#include "stdio.fh"
7#include "util.fh"
8#include "rtdb.fh"
9#include "geom.fh"
10#include "qmmm.fh"
11#include "qmmm_params.fh"
12#include "mm_utils.fh"
13#include "inp.fh"
14c
15      integer rtdb
16c
17      integer ntot
18      integer i_c,h_c
19      integer i_c0,h_c0
20      integer i_ai,h_ai
21      integer i_g,h_g
22      double precision c0,gref
23      integer ia,i,k
24      double precision e0, ep, em, grad, hess, step
25      double precision step0
26      integer i1,k1
27c
28      character*32 pname
29      character*255 tag
30      character*30 region(3)
31      integer nregion
32c
33      pname = "qmmm_check_forces"
34c
35      if (ga_nodeid() .eq. 0) then
36         write(6,*)
37         write(6,*) '@ Checking forces'
38         write(6,*)
39         call util_flush(6)
40      end if
41c
42c
43c     region definitions
44c     ------------------
45      tag ="qmmm:region"
46      if (.not.rtdb_get(rtdb,tag(1:inp_strlen(tag))//"_n",
47     >                 mt_int,1,nregion))
48     >      call errquit(pname//tag,0,RTDB_ERR)
49      if(nregion.gt.3)
50     >      call errquit(pname//"too many regions",0,0)
51      if (.not.rtdb_cget(rtdb,tag,nregion,region))
52     >      call errquit(pname//tag,0,RTDB_ERR)
53c
54c     define set of active atoms
55c     --------------------------
56      call qmmm_cons_reset()
57      call qmmm_cons_set("fix","solute")
58      call qmmm_cons_set("fix","solvent")
59      call qmmm_cons_set("free",region(1))
60c
61
62c
63c     get total number of atoms
64c     -------------------------
65c      call qmmm_cons_get_nacts(ntot)
66      ntot = 1
67      if (ga_nodeid() .eq. 0) then
68         write(6,*)
69         write(6,*) '@ number of atoms',ntot
70         write(6,*)
71         call util_flush(6)
72      end if
73c
74      if(.not.ma_alloc_get(MT_INT, ntot, 'qmmm ref index array',
75     &      h_ai,i_ai) ) call errquit(
76     &      'qmmm_data_alloc: unable to allocate heap space',
77     &      ntot, MA_ERR)
78      call ifill(ntot,-1,int_mb(i_ai),1)
79
80      if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref grad array',
81     &      h_g, i_g) ) call errquit(
82     &      'qmmm_data_alloc: unable to allocate heap space',
83     &      3*ntot, MA_ERR)
84      call dfill(3*ntot,0,dbl_mb(i_g),1)
85
86      if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref coord array',
87     &      h_c,i_c) ) call errquit(
88     &      'qmmm_data_alloc: unable to allocate heap space',
89     &      ntot, MA_ERR)
90      call dfill(3*ntot,0,dbl_mb(i_c),1)
91
92c      call qmmm_cons_get_i_acts(int_mb(i_ai))
93      int_mb(i_ai) = 10
94      call mm_get_solute_coord_raw(ntot,
95     >                       int_mb(i_ai),
96     >                       dbl_mb(i_c))
97
98      call md_sp_qmmm()
99      call prp_print()
100        if (.not. rtdb_get(rtdb,'md:energy', mt_dbl, 1, e0))
101     $     call errquit('driver: could not get energy',0, RTDB_ERR)
102        call mm_get_solute_force_raw(ntot,
103     &                    int_mb(i_ai),
104     &                    dbl_mb(i_g))
105
106      if(.not.rtdb_get(rtdb,'qmmm:step',mt_dbl,
107     >                 1,step0))
108     >   step0=0.01d0
109
110       call mm_print_system()
111
112c
113c     loop over all atoms
114c     -------------------
115      if (ga_nodeid() .eq. 0)
116     >   write(6,8) "comp","anal-g","num-g","error","step","ep","em"
117 8       format(1x,'checkgrad ',A5,4A16,2x,2A12)
118
119      i=0
120      do ia=1,ntot
121        do k=1,3
122          i=i+1
123          c0=dbl_mb(i_c+i-1)
124          gref=dbl_mb(i_g+i-1)
125
126          step = step0
12710        continue
128          dbl_mb(i_c+i-1)=c0+step
129          call mm_set_solute_coord_raw(ntot,
130     >                     int_mb(i_ai),
131     >                     dbl_mb(i_c))
132
133          call md_sp_qmmm()
134          if (.not. rtdb_get(rtdb,'md:energy', mt_dbl, 1, ep))
135     $     call errquit('driver: could not get energy',0, RTDB_ERR)
136          call mm_print_system()
137          call prp_print()
138          write(*,*) "c,c0,ep",dbl_mb(i_c+i-1),c0,ep
139
140c          if (abs(ep-e0) .lt. 1e-6) then
141c             write(6,*) ' Increasing the step ', ep,e0,ep-e0, step
142c             step = step*10.0d0
143c             goto 10
144c          else if (abs(ep-e0) .gt. 1e-2) then
145c             write(6,*) ' Decreasing the step ', ep,e0,ep-e0, step
146c             step = step/3.0
147c             if(step.le.0.00000001d0) goto 12
148c             goto 10
149c          end if
15012        continue
151
152          dbl_mb(i_c+i-1)=c0-step
153          call mm_set_solute_coord_raw(ntot,
154     >                     int_mb(i_ai),
155     >                     dbl_mb(i_c))
156
157          call md_sp_qmmm()
158          if (.not. rtdb_get(rtdb,'md:energy', mt_dbl, 1, em))
159     $     call errquit('driver: could not get energy',0, RTDB_ERR)
160
161          write(*,*) "c,c0,em",dbl_mb(i_c+i-1),c0,em
162
163          grad = (ep - em) / (2.0d0*step)
164c
165          if (ga_nodeid() .eq. 0)
166     >      write(6,7) int_mb(i_ai+ia-1),
167     >         gref, grad,abs(gref-grad), step, ep,em
168 7        format(1x,'checkgrad ',i5,2f16.8,e16.8,f16.8,2x,2f16.8)
169          dbl_mb(i_c+i-1)=c0
170        end do
171      end do
172
173
174      end
175
176      subroutine qmmm_check_forces(rtdb)
177      implicit none
178#include "errquit.fh"
179#include "global.fh"
180#include "mafdecls.fh"
181#include "stdio.fh"
182#include "util.fh"
183#include "rtdb.fh"
184#include "geom.fh"
185#include "qmmm.fh"
186#include "qmmm_params.fh"
187#include "mm_utils.fh"
188#include "inp.fh"
189c
190      integer rtdb
191c
192      integer ntot
193      integer i_c,h_c
194      integer i_c0,h_c0
195      integer i_ai,h_ai
196      integer i_g,h_g
197      double precision c0,gref
198      integer ia,i,k
199      double precision e0, ep, em, grad, hess, step
200      double precision step0
201      integer i1,k1
202c
203      logical qmmm_energy_gradient
204      logical task_qmmm_gradient
205      external qmmm_energy_gradient
206      external task_qmmm_gradient
207      logical task_qmmm_energy
208      external task_qmmm_energy
209      character*32 pname
210      character*255 tag
211      character*30 region(3)
212      integer nregion
213c
214      pname = "qmmm_check_forces"
215c
216      if (ga_nodeid() .eq. 0) then
217         write(6,*)
218         write(6,*) '@ Checking forces'
219         write(6,*)
220         call util_flush(6)
221      end if
222c
223c
224c     region definitions
225c     ------------------
226      tag ="qmmm:region"
227      if (.not.rtdb_get(rtdb,tag(1:inp_strlen(tag))//"_n",
228     >                 mt_int,1,nregion))
229     >      call errquit(pname//tag,0,RTDB_ERR)
230      if(nregion.gt.3)
231     >      call errquit(pname//"too many regions",0,0)
232      if (.not.rtdb_cget(rtdb,tag,nregion,region))
233     >      call errquit(pname//tag,0,RTDB_ERR)
234c
235c     define set of active atoms
236c     --------------------------
237      call qmmm_cons_reset()
238      call qmmm_cons_set("fix","solute")
239      call qmmm_cons_set("fix","solvent")
240      call qmmm_cons_set("free",region(1))
241c
242
243      call mm_print_system()
244c
245c     get total number of atoms
246c     -------------------------
247      call qmmm_cons_get_nacts(ntot)
248      if (ga_nodeid() .eq. 0) then
249         write(6,*)
250         write(6,*) '@ number of atoms',ntot
251         write(6,*)
252         call util_flush(6)
253      end if
254c
255      if(.not.ma_alloc_get(MT_INT, ntot, 'qmmm ref index array',
256     &      h_ai,i_ai) ) call errquit(
257     &      'qmmm_data_alloc: unable to allocate heap space',
258     &      ntot, MA_ERR)
259      call ifill(ntot,-1,int_mb(i_ai),1)
260
261      if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref grad array',
262     &      h_g, i_g) ) call errquit(
263     &      'qmmm_data_alloc: unable to allocate heap space',
264     &      3*ntot, MA_ERR)
265      call dfill(3*ntot,0,dbl_mb(i_g),1)
266
267      if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref coord array',
268     &      h_c,i_c) ) call errquit(
269     &      'qmmm_data_alloc: unable to allocate heap space',
270     &      ntot, MA_ERR)
271      call dfill(3*ntot,0,dbl_mb(i_c),1)
272
273      call qmmm_cons_get_i_acts(int_mb(i_ai))
274      call mm_get_solute_coord(ntot,
275     >                       int_mb(i_ai),
276     >                       dbl_mb(i_c))
277
278      call md_sp_qmmm()
279      if (.not. qmmm_energy_gradient(rtdb,.true.))
280     $    call errquit(pname//'qmmm_energy_gradient failed'
281     $                  ,0, GEOM_ERR)
282        call qmmm_energy_rtdb_push(rtdb)
283        if (.not. rtdb_get(rtdb,'qmmm:energy', mt_dbl, 1, e0))
284     $     call errquit('driver: could not get energy',0, RTDB_ERR)
285        call qmmm_print_energy(rtdb)
286        call mm_get_solute_force(ntot,
287     &                    int_mb(i_ai),
288     &                    dbl_mb(i_g))
289
290      if(.not.rtdb_get(rtdb,'qmmm:step',mt_dbl,
291     >                 1,step0))
292     >   step0=0.01d0
293
294c
295c     loop over all atoms
296c     -------------------
297      if (ga_nodeid() .eq. 0)
298     >   write(6,8) "comp","anal-g","num-g","error","step","ep","em"
299 8       format(1x,'checkgrad ',A5,4A16,2x,2A12)
300
301      i=0
302      do ia=1,ntot
303        do k=1,3
304          i=i+1
305          c0=dbl_mb(i_c+i-1)
306          gref=dbl_mb(i_g+i-1)
307
308          step = step0
30910        continue
310          dbl_mb(i_c+i-1)=c0+step
311          call mm_set_solute_coord(ntot,
312     >                     int_mb(i_ai),
313     >                     dbl_mb(i_c))
314
315          call md_sp_qmmm()
316          if (.not. qmmm_energy_gradient(rtdb,.false.))
317     $    call errquit(pname//'qmmm_energy_gradient failed'
318     $                  ,0, GEOM_ERR)
319          call qmmm_energy_rtdb_push(rtdb)
320          if (.not. rtdb_get(rtdb,'qmmm:energy', mt_dbl, 1, ep))
321     $     call errquit('driver: could not get energy',0, RTDB_ERR)
322          call qmmm_print_energy(rtdb)
323
324c          if (abs(ep-e0) .lt. 1e-6) then
325c             write(6,*) ' Increasing the step ', ep,e0,ep-e0, step
326c             step = step*10.0d0
327c             goto 10
328c          else if (abs(ep-e0) .gt. 1e-2) then
329c             write(6,*) ' Decreasing the step ', ep,e0,ep-e0, step
330c             step = step/3.0
331c             if(step.le.0.00000001d0) goto 12
332c             goto 10
333c          end if
33412        continue
335
336          dbl_mb(i_c+i-1)=c0-step
337          call mm_set_solute_coord(ntot,
338     >                     int_mb(i_ai),
339     >                     dbl_mb(i_c))
340
341          call md_sp_qmmm()
342          if (.not. qmmm_energy_gradient(rtdb,.false.))
343     $    call errquit(pname//'qmmm_energy_gradient failed'
344     $                  ,0, GEOM_ERR)
345          call qmmm_energy_rtdb_push(rtdb)
346          if (.not. rtdb_get(rtdb,'qmmm:energy', mt_dbl, 1, em))
347     $     call errquit('driver: could not get energy',0, RTDB_ERR)
348          call qmmm_print_energy(rtdb)
349
350
351          grad = (ep - em) / (2.0d0*step)
352c
353          if (ga_nodeid() .eq. 0)
354     >      write(6,7) int_mb(i_ai+ia-1),
355     >         gref, grad,abs(gref-grad), step, ep,em
356 7        format(1x,'checkgrad ',i5,2f16.8,e16.8,f16.8,2x,2f16.8)
357          dbl_mb(i_c+i-1)=c0
358        end do
359      end do
360
361
362      end
363
364      subroutine qmmm_check_forces0(rtdb)
365      implicit none
366#include "errquit.fh"
367#include "global.fh"
368#include "mafdecls.fh"
369#include "stdio.fh"
370#include "util.fh"
371#include "rtdb.fh"
372#include "geom.fh"
373#include "qmmm.fh"
374#include "qmmm_params.fh"
375#include "mm_utils.fh"
376c
377      integer rtdb
378c
379      integer ntot
380      integer i_c,h_c
381      integer i_c0,h_c0
382      integer i_ai,h_ai
383      integer i_g,h_g
384      double precision c0,gref
385      integer ia,i,k
386      double precision e0, ep, em, grad, hess, step
387      double precision step0
388      integer i1,k1
389c
390      logical qmmm_energy_gradient
391      external qmmm_energy_gradient
392      logical task_qmmm_energy
393      external task_qmmm_energy
394      character*32 pname
395
396      pname = "qmmm_check_forces"
397c
398       call md_sp()
399       if (.not. rtdb_get(rtdb,'md:energy',mt_dbl,1,e0))
400     $      call errquit('qmmm: failed put energy', 0, RTDB_ERR)
401          e0 = e0/cau2kj
402c
403      if (ga_nodeid() .eq. 0) then
404         write(6,*)
405         write(6,*) '@ Checking forces'
406         write(6,*)
407         call util_flush(6)
408      end if
409c
410c     get total number of atoms
411c     -------------------------
412      call mm_get_solute_tot_na_gen(ntot,mm_link)
413           if (ga_nodeid() .eq. 0) then
414         write(6,*)
415         write(6,*) '@ number of atoms',ntot
416         write(6,*)
417         call util_flush(6)
418      end if
419
420c
421      if(.not.ma_alloc_get(MT_INT, ntot, 'qmmm ref index array',
422     &      h_ai,i_ai) ) call errquit(
423     &      'qmmm_data_alloc: unable to allocate heap space',
424     &      ntot, MA_ERR)
425      call ifill(ntot,-1,int_mb(i_ai),1)
426
427      if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref grad array',
428     &      h_g, i_g) ) call errquit(
429     &      'qmmm_data_alloc: unable to allocate heap space',
430     &      3*ntot, MA_ERR)
431      call dfill(3*ntot,0,dbl_mb(i_g),1)
432
433      if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref coord array',
434     &      h_c,i_c) ) call errquit(
435     &      'qmmm_data_alloc: unable to allocate heap space',
436     &      ntot, MA_ERR)
437      call dfill(3*ntot,0,dbl_mb(i_c),1)
438
439      call mm_get_solute_ind_gen(ntot,
440     >                         mm_link,
441     >                         int_mb(i_ai))
442
443      call mm_get_solute_coord_gen(ntot,
444     >                         mm_link,
445     >                         int_mb(i_ai),
446     >                         dbl_mb(i_c))
447
448      call mm_get_solute_force_gen(ntot,
449     >                         mm_link,
450     >                         int_mb(i_ai),
451     >                         dbl_mb(i_g))
452
453      if(.not.rtdb_get(rtdb,'qmmm:step',mt_dbl,
454     >                 1,step0))
455     >   step0=0.01d0
456
457c      do i=1,3*ntot
458c        dbl_mb(i_g+i-1)=dbl_mb(i_g+i-1)*cau2kj/cau2nm
459c      end do
460c
461c     loop over all atoms
462c     -------------------
463      if (ga_nodeid() .eq. 0)
464     >   write(6,8) "comp","anal-g","num-g","error","step","ep","em"
465 8       format(1x,'checkgrad ',A5,4A16,2x,2A12)
466
467      i=0
468      do ia=1,ntot
469        do k=1,3
470          i=i+1
471          c0=dbl_mb(i_c+i-1)
472          gref=dbl_mb(i_g+i-1)
473
474          step = step0
47510        continue
476          dbl_mb(i_c+i-1)=c0+step
477          call mm_set_solute_coord_gen(ntot,
478     >                       mm_link,
479     >                       int_mb(i_ai),
480     >                       dbl_mb(i_c))
481          call md_sp()
482          if (.not. rtdb_get(rtdb,'md:energy',mt_dbl,1,ep))
483     $      call errquit('qmmm: failed put energy', 0, RTDB_ERR)
484          ep = ep/cau2kj
485
486          if (abs(ep-e0) .lt. 1e-6) then
487             write(6,*) ' Increasing the step ', ep,e0,ep-e0, step
488             step = step*10.0d0
489             goto 10
490          else if (abs(ep-e0) .gt. 1e-2) then
491             write(6,*) ' Decreasing the step ', ep,e0,ep-e0, step
492             step = step/3.0
493             if(step.le.0.00000001d0) goto 12
494             goto 10
495          end if
49612        continue
497
498          dbl_mb(i_c+i-1)=c0-step
499          call mm_set_solute_coord_gen(ntot,
500     >                       mm_link,
501     >                       int_mb(i_ai),
502     >                       dbl_mb(i_c))
503c
504            call md_sp()
505          if (.not. rtdb_get(rtdb,'md:energy',mt_dbl,1,em))
506     $      call errquit('qmmm: failed put energy', 0, RTDB_ERR)
507          em = em/cau2kj
508c
509          grad = (ep - em) / (2.0d0*step)
510c
511          if (ga_nodeid() .eq. 0)
512     >      write(6,7) int_mb(i_ai+ia-1),
513     >         gref, grad,abs(gref-grad), step, ep,em
514 7        format(1x,'checkgrad ',i5,2f16.8,e16.8,f16.8,2x,2f16.8)
515          dbl_mb(i_c+i-1)=c0
516        end do
517      end do
518
519
520      end
521
522c $Id$
523