1      subroutine qmmm_forces(inrtdb,mwm,nwm,mwa,nwa,iwz,xw,fw,
2     + msa,nsa,isat,isdt,isq,isz,xs,fs,energy)
3c
4c $Id$
5c
6      implicit none
7#include "errquit.fh"
8c
9#include "rtdb.fh"
10#include "geom.fh"
11#include "util.fh"
12#include "global.fh"
13#include "mafdecls.fh"
14#include "qmmm.fh"
15#include "qmmm_params.fh"
16c
17c
18      integer inrtdb
19      integer mwm,nwm,mwa,nwa,msa,nsa
20      integer iwz(mwm),isat(msa),isdt(msa),isq(msa),isz(msa)
21      real*8 energy,xw(mwm,3,mwa),fw(mwm,3,mwa),xs(msa,3),fs(msa,3)
22c
23      character*32 interface
24      character*32 optimization
25      character*255 theory
26      character*30 operation
27      logical qmmm_init
28      character*32 pname
29      double precision eatoms
30      logical status
31      logical do_grad
32
33      logical qmmm_energy_gradient
34      external qmmm_energy_gradient
35
36      pname = 'qmmm_main'
37      if(qmmm_print_debug())
38     >  write(*,*) "in "//pname
39c
40c     make sure qmmm is initialized
41c     -----------------------------
42      if (.not. rtdb_get(inrtdb,'qmmm:init',mt_log,1,qmmm_init))
43     $     qmmm_init=.false.
44      if(.not.qmmm_init) return
45c
46c     nothing should be calculated if
47c     qm module (e.g. driver is in charge)
48c     -----------------------------------
49      interface = qmmm_get_interface()
50      if(interface.eq.'qm') then
51      if(qmmm_print_debug())
52     > write(*,*) "exiting out of qmmm_forces since interface is set qm"
53        energy = 0
54        return
55      end if
56
57      if (.not. rtdb_cget(inrtdb, 'qmmm:operation', 1, operation))
58     $     operation = ' '
59
60      if(operation.eq.'energy') then
61        do_grad=.false.
62      else
63        do_grad=.true.
64      end if
65
66      call  qmmm_bq_data_reload()
67      if(.not.qmmm_energy_gradient(inrtdb,do_grad))
68     $   call errquit(pname//'failed qmmm_energy_gradient',0,0)
69
70      if(do_grad) call qmmm_cons_fixed()
71
72      if (.not.rtdb_get(inrtdb,'qmmm:uqmatm',mt_dbl,1,eatoms))
73     $     call errquit('qmmm: failed getting  ref energy',0,RTDB_ERR)
74
75      if (.not. rtdb_get(inrtdb,'qmmm:qm_energy',mt_dbl,1,energy))
76     $     call errquit('qmmm: failed getting qm energy',0,RTDB_ERR)
77
78      energy = energy-eatoms
79      energy=energy*cau2kj
80
81      return
82      end
83
84