1      logical function drv_lst(rtdb)
2*
3* $Id$
4*
5      implicit none
6#include "mafdecls.fh"
7#include "rtdb.fh"
8#include "global.fh"
9#include "stdio.fh"
10#include "util.fh"
11      integer  rtdb
12      logical  hnd_lstx
13      external hnd_lstx
14      integer  ir, iw
15      logical  some, dbug
16      common/hnd_iofile/ir,iw
17c
18      call util_print_push
19      call util_print_rtdb_load(rtdb, 'driver')
20      call ecce_print_module_entry('driver')
21      if (util_print('brdcst', print_never)) call setdbg(1)
22c
23      dbug=.false.
24      some=.false.
25      some=some.or.dbug
26      ir=LuIn
27      iw=LuOut
28      if(some) then
29         write(iw,*) ' drv_lst: calling hnd_lstx . '
30      endif
31c
32c     ----- lst path -----
33c
34      drv_lst = hnd_lstx(rtdb)
35c
36      if (drv_lst) then
37         call ecce_print_module_exit('driver', 'ok')
38      else
39         call ecce_print_module_exit('driver', 'failed')
40      endif
41c
42      if (util_print('brdcst', print_never)) call setdbg(0)
43      call util_print_pop
44c
45      end
46      logical function hnd_lstx(rtdb)
47      implicit double precision (a-h,o-z)
48#include "errquit.fh"
49#include "rtdb.fh"
50#include "msgtypesf.h"
51#include "global.fh"
52#include "mafdecls.fh"
53#include "stdio.fh"
54#include "geom.fh"
55#include "util.fh"
56      logical       geom_zmt_get_nzvar
57      logical       geom_zmt_get_nizmat
58      external      geom_zmt_get_nizmat
59      logical       geom_zmt_get_nzfrz
60      external      geom_zmt_get_nzfrz
61      logical       geom_zmt_get_izfrz
62      external      geom_zmt_get_izfrz
63      logical       geom_lst_get_coord
64      external      geom_lst_get_coord
65      integer       rtdb
66      integer       geom
67      logical       do_gradient
68      logical       some
69      logical       out
70      logical       dbug
71      logical       status
72      logical       done
73      logical       rstart
74      logical       zcoord
75      character*16  atmlab
76      parameter (mxatom=500)
77      parameter (mxcart=3*mxatom)
78      parameter (mxzmat=1500)
79      parameter (mxcoor=1500)
80c                mxcoor=max(mxcart,mxzmat)
81      common/hnd_iofile/ir,iw
82      common/hnd_optmiz/x0(mxcoor),x(mxcoor),dx(mxcoor),
83     1                  g0(mxcoor),g(mxcoor),ds(mxcoor),
84     2                  func,func0,gmax,gmax0,curv,alpha,gnorm
85      common/hnd_optfrz/nzfrz,izfrz(mxcoor),iatfrz(mxatom)
86      common/hnd_optvar/zcoord,ncoord,mcoord
87      common/hnd_optfun/e,eg(mxcart)
88      common/hnd_lstrun/nptlst,iptlst
89      common/hnd_lsttim/energy_time,gradient_time
90      common/hnd_zmtpar/nzmat,nzvar,nvar
91      common/hnd_molxyz/c(mxcart),zan(mxatom),nat
92      common/hnd_mollab/atmlab(mxatom)
93      data zero   /0.0d+00/
94      data pt5    /0.5d+00/
95      data one    /1.0d+00/
96c
97      hnd_lstx = .false.
98c
99      do_gradient = .false.
100c
101      dbug=.false.
102      out =.false.
103      out =out.or.dbug
104      some=.false.
105      some=some.or.out
106      if (ga_nodeid().eq.0.or.some)  then
107         write(iw,9999)
108      endif
109c
110c     ----- get going ... -----
111c
112      if (.not. geom_create(geom, 'geometry'))
113     &       call errquit('hnd_lst: geom_create?', 911, GEOM_ERR)
114      if (.not. geom_rtdb_load(rtdb, geom, 'geometry'))
115     &       call errquit('hnd_lst: no geometry ', 911, RTDB_ERR)
116      if (.not. geom_cart_get(geom, ncent, atmlab, c, zan))
117     &       call errquit('hnd_lst: geom_get ', 911, GEOM_ERR)
118       nat = ncent
119      if(.not.geom_zmt_get_nizmat(geom,nzmat))
120     &       call errquit('geom_input: geom_zmt_get_nizmat failed',0,
121     &       GEOM_ERR)
122      if(.not.geom_zmt_get_nzvar(geom,nzvar))
123     &     call errquit('geom_input: geom_zmt_get_nzvar failed',0,
124     &       GEOM_ERR)
125c
126c     ----- scan -lst- coordinates file -----
127c
128      nptlst=1
129   10 continue
130         status=geom_lst_get_coord(x(1),x(1+nat),x(1+nat),
131     &                             nat,nptlst)
132         if(status) then
133            nptlst=nptlst+1
134            go to 10
135         else
136            nptlst=nptlst-1
137         endif
138c
139      if (ga_nodeid().eq.0.or.dbug)  then
140         write(iw,9996) nptlst
141      endif
142      if(nptlst.eq.0) then
143         write(iw,*) 'no -lst- points found. return and continue ..'
144         return
145      endif
146c
147      nzfrz=0
148      do i=1,mxcoor
149         izfrz(i)=0
150      enddo
151      do i=1,mxatom
152         iatfrz(i)=0
153      enddo
154c
155c     ----- frozen coordinates ? if so .... -----
156c
157      nzfrz=0
158      if (geom_zmt_get_nzfrz(geom,nzfrz)) then
159         if(.not.geom_zmt_get_izfrz(geom,izfrz,nzfrz))
160     $        call errquit('geom_input: geom_zmt_get_izfrz  failed',0,
161     &       GEOM_ERR)
162      endif
163      if (ga_nodeid().eq.0.or.dbug)  then
164         if(nzfrz.gt.0) then
165            write(iw,9993) nzfrz
166            write(iw,9992) (izfrz(i),i=1,nzfrz)
167         else
168            write(iw,9991)
169         endif
170      endif
171c
172c     ----- get ready to start now -----
173c
174c     if internal coordinates ...
175c           nzvar  = # of (redundant) internal coordinates
176c           ncoord = # of independent internal coordinates
177c           mcoord = max( 3*nat , nzvar )
178c     if cartesian coordinates ...
179c           ncoord = 3*nat
180c           mcoord = 3*nat
181c
182      rstart     =.false.
183c
184      zcoord=.false.
185      zcoord=nzmat.gt.0
186      if(dbug)  then
187         write(iw,*) 'nzmat,zcoord = ',nzmat,zcoord
188         write(iw,*) 'nzfrz        = ',nzfrz
189         if(nzfrz.gt.0) then
190            write(iw,*) 'izfrz = ',(izfrz(i),i=1,nzfrz)
191         endif
192      endif
193      if(zcoord) then
194         call geom_bandbi(geom)
195         ncart =3*nat
196         ncoord=nvar
197         mcoord=max(ncart,nzvar)
198      else
199         nzmat =0
200         nzvar =0
201         nvar  =3*nat
202         ncart =3*nat
203         ncoord=ncart
204         mcoord=ncart
205      endif
206      if(dbug)  then
207         write(iw,9997) nzvar,nvar,ncart,ncoord,mcoord
208      endif
209c
210c     ----- set up ----
211c
212      e     =zero
213      func  =zero
214      func0 =zero
215      gmax  =zero
216      gmax0 =zero
217      curv  =zero
218      do i=1,mxcoor
219         dx(i)=zero
220         ds(i)=zero
221         eg(i)=zero
222         g0(i)=zero
223         x0(i)=zero
224          x(i)=zero
225          g(i)=zero
226      enddo
227c
228c     ----- first point -----
229c
230      iptlst=1
231      do i=1,ncart
232         x(i)=c(i)
233      enddo
234c
235c     ----- energy -----
236c
237         time_start=util_wallsec()
238      call hnd_lst_energy(rtdb,geom)
239         time_end  =util_wallsec()
240         energy_time=time_end-time_start
241         if(out) then
242            write(iw,*) 'energy_time = ',energy_time
243         endif
244c
245      func=e
246      if(do_gradient) then
247c
248c     ----- gradient -----
249c
250            time_start=util_wallsec()
251         call hnd_opt_gradient(rtdb,geom)
252            time_end  =util_wallsec()
253            gradient_time=time_end-time_start
254            if(out) then
255               write(iw,*) 'gradient_time = ',gradient_time
256            endif
257c
258         do i=1,ncart
259            g(i)=eg(i)
260         enddo
261         if(zcoord) then
262            call hnd_opt_tfgx(rtdb,geom)
263         endif
264      endif
265c
266      call hnd_lst_print(rtdb)
267      do i=1,ncart
268         x0(i)=x(i)
269      enddo
270      do i=1,ncoord
271         g0(i)=g(i)
272      enddo
273c
274      if (ga_nodeid().eq.0.and.dbug)  then
275         if (.not. geom_print(geom)) call errquit
276     $        ('hnd_lst_drv: geom_print?',0, GEOM_ERR)
277         if (util_print('bonds',print_default)) then
278           if (.not.geom_print_distances(geom)) call errquit(
279     &         'hnd_lst_drv: geom_print_distances failed',911,
280     &       GEOM_ERR)
281         endif
282         if (util_print('angles',print_default)) then
283           if (.not.geom_print_angles(geom)) call errquit(
284     &         'hnd_lst_drv: geom_print_angles failed',911, GEOM_ERR)
285         endif
286      endif
287      call ga_sync()
288c
289c     ----- lst path -----
290c
291  100 continue
292      iptlst=iptlst+1
293c
294      call hnd_lst_path(rtdb,geom,rstart)
295c
296      if(zcoord) then
297         call hnd_opt_tfgx(rtdb,geom)
298      endif
299      call hnd_lst_print(rtdb)
300c
301      if (ga_nodeid().eq.0.and.dbug)  then
302         if (.not. geom_print(geom)) call errquit
303     $        ('hnd_lst_drv: geom_print?',0, GEOM_ERR)
304         if (util_print('bonds',print_default)) then
305           if (.not.geom_print_distances(geom)) call errquit(
306     &         'hnd_lst_drv: geom_print_distances failed',911, GEOM_ERR)
307         endif
308         if (util_print('angles',print_default)) then
309           if (.not.geom_print_angles(geom)) call errquit(
310     &         'hnd_lst_drv: geom_print_angles failed',911,
311     &       GEOM_ERR)
312         endif
313      endif
314      call ga_sync()
315c
316c     ----- go back to next -lst- point -----
317c
318      done=iptlst.eq.nptlst
319      if(.not.done) then
320         go to 100
321      else
322         if (ga_nodeid().eq.0.or.dbug)  then
323            write(iw,9998)
324         endif
325         call ga_sync()
326         if (.not.geom_destroy(geom))
327     &      call errquit('hnd_lst: geom_destroy?',911, GEOM_ERR)
328         hnd_lstx=.true.
329      endif
330c
331 9999 format(/,10x,8(1h-),
332     1       /,10x,'LST Path',
333     2       /,10x,8(1h-))
334 9998 format(/,10x,18(1h-),
335     1       /,10x,'LST Path completed',
336     2       /,10x,18(1h-))
337 9997 format(' in hnd_lstx, nzvar,nvar,ncart,ncoord,mcoord = ',5i5)
338 9996 format(' number of points on -lst- path (nptlst) = ',i4)
339 9993 format(' number of frozen internal coordinates = ',i4)
340 9992 format('           frozen internal coordinates = ',12i4)
341 9991 format('        no frozen internal coordinates.')
342      end
343      SUBROUTINE HND_LST_PRINT(RTDB)
344      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
345#include "mafdecls.fh"
346#include "rtdb.fh"
347#include "global.fh"
348      INTEGER   RTDB
349      PARAMETER (MXATOM=500)
350      PARAMETER (MXCART=3*MXATOM)
351      PARAMETER (MXZMAT=1500)
352      PARAMETER (MXCOOR=1500)
353      LOGICAL      DBUG
354      CHARACTER*16 ATMNAM
355      COMMON/HND_IOFILE/IR,IW
356      COMMON/HND_MOLLAB/ATMNAM(MXATOM)
357      COMMON/HND_MOLNUC/NUC(MXATOM)
358      COMMON/HND_MOLXYZ/C(3,MXATOM),ZAN(MXATOM),NAT
359      COMMON/HND_OPTMIZ/X0(MXCOOR),X(MXCOOR),DX(MXCOOR),
360     1                  G0(MXCOOR),G(MXCOOR),DS(MXCOOR),
361     2                  FUNC,FUNC0,GMAX,GMAX0,CURV,ALPHA,GNORM
362      COMMON/HND_LSTRUN/NPTLST,IPTLST
363      COMMON/HND_OPTFRZ/NZFRZ,IZFRZ(MXCOOR),IATFRZ(MXATOM)
364      COMMON/HND_OPTFUN/E,EG(MXCART)
365C
366      DBUG=.FALSE.
367C
368C     ----- PRINT OPTIMIZATION SUMMARY -----
369C
370      IF( GA_NODEID().EQ.0.OR.DBUG) THEN
371         WRITE(IW,9999) NPTLST,IPTLST,FUNC
372      ENDIF
373      RETURN
374 9999 FORMAT(1H1,/,1X,
375     1       'nptlst  iptlst      func             ',1X,I5,I8,F17.8)
376      END
377      SUBROUTINE HND_LST_PATH(RTDB,GEOM,RSTART)
378      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
379#include "mafdecls.fh"
380#include "tcgmsg.fh"
381#include "msgtypesf.h"
382#include "global.fh"
383#include "geom.fh"
384#include "rtdb.fh"
385C
386C     ----- ONE DIMENSIONAL SEARCH.  -----
387C
388      PARAMETER (MXATOM=500)
389      PARAMETER (MXCART=3*MXATOM)
390      PARAMETER (MXZMAT=1500)
391      PARAMETER (MXCOOR=1500)
392C
393      INTEGER   RTDB
394      INTEGER   GEOM
395      LOGICAL   STATUS
396      LOGICAL   UTIL_TEST_TIME_REMAINING
397      EXTERNAL  UTIL_TEST_TIME_REMAINING
398      LOGICAL   GEOM_LST_GET_COORD
399      EXTERNAL  GEOM_LST_GET_COORD
400      CHARACTER*16 TAGS_NW
401      DIMENSION COORDS_NW(MXCART)
402      DIMENSION CHARGE_NW(MXATOM)
403      DIMENSION   TAGS_NW(MXATOM)
404C
405      LOGICAL DO_GRADIENT
406      LOGICAL RSTART
407      LOGICAL DBUG
408      LOGICAL OUT
409      COMMON/HND_IOFILE/IR,IW
410      COMMON/HND_MOLXYZ/C(MXCART),ZAN(MXATOM),NAT
411      COMMON/HND_OPTMIZ/X0(MXCOOR),X(MXCOOR),DX(MXCOOR),
412     1                  G0(MXCOOR),G(MXCOOR),DS(MXCOOR),
413     2                  FUNC,FUNC0,GMAX,GMAX0,DELSQ,ALPHA,GNORM
414      COMMON/HND_LSTRUN/NPTLST,IPTLST
415      COMMON/HND_OPTFRZ/NZFRZ,IZFRZ(MXCOOR),IATFRZ(MXATOM)
416      COMMON/HND_OPTTIM/ENERGY_TIME,GRADIENT_TIME
417      COMMON/HND_OPTFUN/E,EG(MXCART)
418      COMMON/HND_ZMTPAR/NZMAT,NZVAR,NVAR
419      DATA ZERO,TWO,THREE       /0.0D+00,2.0D+00,3.0D+00/
420      DATA ONEPT5               /1.5D+00/
421C
422      DO_GRADIENT=.FALSE.
423C
424      DBUG=.FALSE.
425      OUT =.FALSE.
426      OUT =OUT.OR.DBUG
427      IF(OUT) THEN
428         WRITE(IW,9994)
429      ENDIF
430C
431      NCART =3*NAT
432C
433C     ----- NORMAL SETUP -----
434C
435      FUNC0 = FUNC
436      DO I = 1,NCART
437         X0(I) = X(I)
438         G0(I) =EG(I)
439         G (I) =EG(I)
440         X (I) =ZERO
441      ENDDO
442C
443C     ----- NEW -LST- POINT -----
444C
445      STATUS=GEOM_LST_GET_COORD(X(1),X(1+NAT),X(1+NAT*2),
446     1                          NAT,IPTLST)
447      DO IAT = 1,NAT
448         C(1+3*(IAT-1)) = X(IAT)
449         C(2+3*(IAT-1)) = X(IAT+NAT)
450         C(3+3*(IAT-1)) = X(IAT+NAT*2)
451      ENDDO
452      DO I = 1,NCART
453         X(I)=ZERO
454         X(I)=C(I)
455      ENDDO
456C
457C     ----- WRITE TO -NWCHEM- -----
458C
459      STATUS=GEOM_CART_GET(GEOM,NAT_NW,TAGS_NW,COORDS_NW,
460     1                                         CHARGE_NW)
461C
462      IF(DBUG) THEN
463         WRITE(IW,9997)
464         WRITE(IW,9998) (       X0(I),I=1,NCART)
465         WRITE(IW,9996)
466         WRITE(IW,9998) (        X(I),I=1,NCART)
467         WRITE(IW,9999)
468         WRITE(IW,9998) (COORDS_NW(I),I=1,NCART)
469      ENDIF
470      DO I = 1,NCART
471         COORDS_NW(I)=X(I)
472      ENDDO
473      STATUS=GEOM_CART_SET(GEOM,NAT_NW,TAGS_NW,COORDS_NW,
474     1                                         CHARGE_NW)
475      STATUS=GEOM_RTDB_STORE(RTDB,GEOM,'geometry')
476C
477      IF(OUT) THEN
478         WRITE(IW,9999)
479         WRITE(IW,9998) (COORDS_NW(I),I=1,NCART)
480      ENDIF
481C
482C
483C     ----- CALL FUNCTION EVALUATION -----
484C
485         STATUS=UTIL_TEST_TIME_REMAINING(RTDB,INT(ENERGY_TIME*ONEPT5))
486         IF(OUT) THEN
487            WRITE(IW,*)
488     1      'ENOUGH TIME REMAINING FOR ENERGY = ? ',STATUS
489         ENDIF
490         IF(.NOT.STATUS) THEN
491            WRITE(IW,*)
492     1      'NOT ENOUGH TIME REMAINING, SHUTTING DOWN ... '
493            RETURN
494         ENDIF
495C
496      CALL HND_LST_ENERGY(RTDB,GEOM)
497      FUNC=E
498C
499C     ----- CALCULATE THE GRADIENT FOR THE FINAL POINT -----
500C
501  900 CONTINUE
502      IF(DO_GRADIENT) THEN
503         STATUS=UTIL_TEST_TIME_REMAINING(RTDB,INT(GRADIENT_TIME*ONEPT5))
504         IF(OUT) THEN
505            WRITE(IW,*)
506     1      'ENOUGH TIME REMAINING FOR GRADIENT = ? ',STATUS
507         ENDIF
508         IF(.NOT.STATUS) THEN
509            WRITE(IW,*)
510     1      'NOT ENOUGH TIME REMAINING, SHUTTING DOWN ... '
511            RETURN
512         ENDIF
513C
514         CALL HND_OPT_GRADIENT(RTDB,GEOM)
515         DO I=1,NCART
516            G(I) =EG(I)
517         ENDDO
518      ENDIF
519C
520C     ----- RETURN AFTER SUCCESSFUL POINT -----
521C
522 1000 CONTINUE
523      IF(DBUG) THEN
524         WRITE(IW,9993)
525      ENDIF
526      RETURN
527 9999 FORMAT(' IN SEARCH_LIN, COORDS = ')
528 9998 FORMAT(F12.7)
529 9997 FORMAT(' IN SEARCH_LIN, X0 = ')
530 9996 FORMAT(' IN SEARCH_LIN, DX = ')
531 9995 FORMAT(' IN SEARCH_LIN,  X = ')
532 9994 FORMAT(' STARTING SEARCH_LIN ')
533 9993 FORMAT('   ENDING SEARCH_LIN ')
534      END
535      subroutine hnd_lst_energy(rtdb,geom)
536      implicit double precision (a-h,o-z)
537#include "errquit.fh"
538#include "mafdecls.fh"
539#include "rtdb.fh"
540#include "geom.fh"
541#include "global.fh"
542      integer  rtdb
543      integer  geom
544      logical  task_energy
545      external task_energy
546      logical  status
547      logical  dbug
548      parameter (mxatom=500)
549      parameter (mxcart=3*mxatom)
550      parameter (mxzmat=1500)
551      parameter (mxcoor=1500)
552      common/hnd_iofile/ir,iw
553      common/hnd_lstrun/nptlst,iptlst
554      common/hnd_optfun/e,eg(mxcart)
555      common/hnd_molxyz/c(3,mxatom),zan(mxatom),nat
556      common/hnd_zmtpar/nzmat,nzvar,nvar
557c
558      logical geom_print_zmatrix
559c
560      dbug=.false.
561c
562      if( ga_nodeid().eq.0.or.dbug) then
563         write(iw,9999) nptlst,iptlst
564         if(dbug) then
565            write(iw,9996)
566            do iat=1,nat
567               write(iw,9995) iat,(c(i,iat),i=1,3)
568            enddo
569            write(iw,9997)
570         endif
571         if(.not.geom_print(geom))
572     1     call errquit('hnd_lst_energy: print error',911, GEOM_ERR)
573      endif
574c
575      if(nzmat.gt.0) then
576         call geom_bandbi(geom)
577         if (.not. geom_print_zmatrix(geom,0d0,.false.))
578     $        call errquit('hnd_lst: print zmat failed',0, GEOM_ERR)
579         write(iw,9997)
580      endif
581c
582      if (task_energy(rtdb)) then
583         status=rtdb_get(rtdb,'task:energy',MT_DBL,1,e)
584         if(dbug) then
585            write(iw,9998) e
586         endif
587      else
588         call errquit('lst path: energy failed', 0, GEOM_ERR)
589      endif
590c
591      return
592 9999 format(1h1,' nptlst',i3,//,'  iptlst',I3)
593 9998 format(' in lst_energy, e = ',f15.10)
594 9997 format(/)
595 9996 format(21x,21(1h-),/,
596     1       21x,'cartesian coordinates',/,
597     2       21x,21(1h-))
598 9995 format(9x,i5,3f15.8)
599      end
600