1      logical function lbfgs_driver(irtdb)
2      implicit none
3#include "mafdecls.fh"
4#include "errquit.fh"
5#include "geom.fh"
6#include "rtdb.fh"
7#include "global.fh"
8      integer irtdb
9      double precision eps,xtol,gtol,stpmin,stpmax
10      integer iprint(2),iflag,icall,n,m,mp,lp,j,iter
11      logical diagco
12      integer master
13c
14c     the driver for LBFGS must always declare LB2 as EXTERNAL
15c
16      external lb2
17      common /lb3/mp,lp,gtol,stpmin,stpmax
18      common /lb3p/master
19
20      integer igeom
21      integer ncent
22      integer i_c,h_c
23      integer i_g,h_g
24      integer i_scr,h_scr
25      integer i_d,h_d
26      integer nscr
27      integer maxiter,nrest
28      double precision energy,stp1
29      logical task_gradient
30      external task_gradient
31      character*32 pname
32
33      pname = "lbfgs_driver"
34
35      master=ga_nodeid()
36      m=5
37      diagco= .false.
38      eps= 1.0d-5
39      xtol= 1.0d-16
40      icall=0
41
42      if(.not.geom_create(igeom,'geometry'))
43     + call errquit(pname//'Failed to create geometry',0, GEOM_ERR)
44
45      if(.not.geom_rtdb_load(irtdb,igeom,'geometry'))
46     + call errquit(pname//'Failed to load geometry',0, GEOM_ERR)
47
48      if(.not. geom_ncent(igeom,ncent) )
49     >    call errquit("qmmm:geom_ncent",0,0)
50
51      if(.not.ma_alloc_get(mt_dbl,3*ncent,'c',h_c,i_c))
52     + call errquit( pname//'Failed to allocate memory for c',
53     + 3*ncent, MA_ERR)
54
55      if(.not.ma_alloc_get(mt_dbl,3*ncent,'g',h_g,i_g))
56     + call errquit( pname//'Failed to allocate memory for c',
57     + 3*ncent, MA_ERR)
58
59      nscr = 3*ncent*(2*m+1)+2*m
60      if(.not.ma_alloc_get(mt_dbl,nscr,'scratch',h_scr,i_scr))
61     + call errquit( pname//'Failed to allocate memory for scr',
62     + 3*ncent, MA_ERR)
63
64      if(.not.ma_alloc_get(mt_dbl,3*ncent,'diag',h_d,i_d))
65     + call errquit( pname//'Failed to allocate memory for scr',
66     + 3*ncent, MA_ERR)
67
68       if(.not.geom_cart_coords_get(igeom,dbl_mb(i_c)))
69     + call errquit('qmmm: Failed to get coord ',0, GEOM_ERR)
70
71      if(.not.rtdb_get(irtdb,'lbfgs:maxiter',mt_int,1,maxiter))
72     +  maxiter = 40
73
74c
75      stp1 = 0.0d0
7610    continue
77      n=100
78      iprint(1)= 1
79      iprint(2)= 0
80      iflag=0
81c
82      lbfgs_driver = .true.
83 20   continue
84      if (.not. task_gradient(irtdb))
85     $   call errquit('driver: task_gradient failed',0, GEOM_ERR)
86      if (.not. rtdb_get(irtdb,'task:energy', mt_dbl, 1, energy))
87     $   call errquit('driver: could not get energy',0, RTDB_ERR)
88      if (.not. rtdb_get(irtdb, 'task:gradient', mt_dbl, 3*ncent,
89     $  dbl_mb(i_g))) call errquit('driver: could not get gradient',0,0)
90
91      write(*,*) "stp1, energy",stp1, energy
92      call lbfgs(3*ncent,m,dbl_mb(i_c),energy,
93     >           dbl_mb(i_g),diagco,
94     >           dbl_mb(i_d),iprint,eps,xtol,
95     >           dbl_mb(i_scr),iflag,iter,stp1)
96
97c
98c     if line search exited due to may function evaluations
99c     try restarting
100      if(iflag.eq.-1) then
101       if(master.eq.0) then
102          write(6,*) "@ too many function evaluations"
103          write(6,*) "@ restarting the optimization"
104          call util_flush(6)
105       end if
106c       iflag=0
107c       call dfill(nscr,0.0d0,dbl_mb(i_scr),1)
108c       stp1=1.0
109       go to 20
110      end if
111      if(iflag.le.0) then
112          lbfgs_driver=.false.
113          go to 50
114      end if
115      icall=icall + 1
116      if(iter.gt.maxiter) go to 50
117      if(.not.geom_cart_coords_set(igeom,dbl_mb(i_c)))
118     + call errquit('qmmm: Failed to get coord ',0, GEOM_ERR)
119      if(.not.geom_rtdb_store(irtdb,igeom,'geometry'))
120     + call errquit(pname//'Failed to load geometry',0, GEOM_ERR)
121      go to 20
122  50  continue
123
124      if(.not.ma_free_heap(h_d))
125     & call errquit(pname//'Failed to deallocate stack c',ncent,
126     &       MA_ERR)
127
128
129      if(.not.ma_free_heap(h_scr))
130     & call errquit(pname//'Failed to deallocate stack c',ncent,
131     &       MA_ERR)
132
133      if(.not.ma_free_heap(h_g))
134     & call errquit(pname//'Failed to deallocate stack c',ncent,
135     &       MA_ERR)
136
137
138      if(.not.ma_free_heap(h_c))
139     & call errquit(pname//'Failed to deallocate stack c',ncent,
140     &       MA_ERR)
141
142       if(.not.geom_destroy(igeom))
143     + call errquit('qmmm: Failed to destroy geometry',0, GEOM_ERR)
144
145      end
146c
147c     ** last line of simple driver (sdrive) **
148
149c $Id$
150