1c $Id: bq_data.F 23019 2012-10-30 00:59:12Z d3y133 $
2
3      function  mmi_init(rtdb)
4      implicit none
5#include "util.fh"
6#include "errquit.fh"
7#include "inp.fh"
8#include "stdio.fh"
9#include "rtdb.fh"
10#include "mafdecls.fh"
11#include "mm_data.fh"
12      logical mmi_init
13      integer rtdb
14c
15      logical ignore
16      logical scnb_adj
17      character*30 operation
18      character*255 crdparmfile
19      character*255 field
20      character*84 tag
21      integer fn
22      character*30 pname
23
24      pname = "mm_init"
25c     write(*,*) pname
26
27      ignore = MA_set_hard_fail(.true.)
28      ignore=MA_set_auto_verify(.true.)
29
30      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
31
32      if(operation.eq."neb") then
33        call mm_neb_init(rtdb)
34      else
35        call mm_coords_init(rtdb)
36        call mm_bndparm_init(rtdb)
37        call mm_bond_coords_init()
38        call mm_vdw_init(rtdb)
39        call mm_vdw_coords_init()
40c       call mm_bond_call_test()
41        call mm_links_init(rtdb)
42
43        scnb_adj = .true.
44        if(.not.rtdb_get(rtdb,'mm:scnb_adjust',mt_log,1,scnb_adj))
45     >    scnb_adj = .false.
46c       adjust scnb if true
47        if(scnb_adj) call mm_vdw_adj_scnb()
48
49        call mm_geom_init(rtdb,"geometry")
50        call mm_bq_init(rtdb)
51      end if
52
53      mmi_init = .true.
54      return
55911   mmi_init = .false.
56      return
57      end
58
59      function  mmi_end(rtdb)
60      implicit none
61#include "rtdb.fh"
62#include "mafdecls.fh"
63#include "global.fh"
64#include "util.fh"
65#include "mm_link_data.fh"
66#include "mm_geom_data.fh"
67      integer rtdb
68
69      logical mmi_end
70
71      logical task_energy
72      external task_energy
73c
74
75      character*30 pname
76      character*30 operation
77      logical ignore
78      logical push_geom
79      logical oprint
80      logical rtdb_mode
81      integer i, k
82      integer master
83      character*16 tag
84      character*50 filename
85      character*255 full_filename
86
87      pname = "mm_end"
88c     write(*,*) pname
89
90      master=0
91
92c     TP: if operation = "optimize" run final single point energy
93c     and  aux_geom = .false.
94      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
95      push_geom = .false.
96      if(operation.eq."optimize") push_geom = .true.
97      if(operation.eq."saddle") push_geom = .true.
98
99      if(push_geom) then
100        aux_geom = .false.
101        call mm_geom_push_active(rtdb)
102        call mm_geom_create_full(rtdb)
103        call mm_link_update_bq_coords(rtdb)
104      end if
105
106      rtdb_mode = rtdb_parallel(.false.)
107
108c     Let master node handle for writing necessary outputs
109      if(ga_nodeid().eq.master) then
110        oprint = .true.
111        if(.not.rtdb_get(rtdb,"qmmm:active_region_xyz",mt_log,1,oprint))
112     >   oprint = .false.
113
114
115        if(oprint) then
116          call util_file_prefix('act.xyz',filename)
117          call util_file_name_noprefix(filename,.false.,
118     >                                 .false.,
119     >                                  full_filename)
120          call mm_write_active_region_xyz(24, full_filename)
121        end if
122
123c       Write out final crdparms file.
124c       If opeation=neb crdparms file for each bead is written
125c       in mm_add_egrad.
126        if(operation.eq."optimize".or.operation.eq."saddle") then
127          call mm_create_crdparm(rtdb)
128        end if
129      end if !! if(ga_nodeid().eq.master)
130
131      call ga_sync()
132
133      rtdb_mode = rtdb_parallel(rtdb_mode)
134
135      ignore = rtdb_delete(rtdb,'bq')
136
137c     Deallocation!!
138      call mm_coords_deallocate()
139      call mm_vdw_deallocate()
140      call mm_vdw_coords_deallocate()
141      call mm_bonded_deallocate()
142      call mm_bond_coords_end()
143      call mm_links_end()
144      call mm_geom_end()
145      call mm_bq_end()
146
147      mmi_end = .true.
148      return
149911   mmi_end = .false.
150      return
151      end
152
153      subroutine mm_test(n,t,c)
154      implicit none
155
156      integer n
157      character*(16) t(n)
158      double precision c(3,n)
159      integer i
160
161      do i=1,n
162         write(6,*) t(i),c(1,i),c(2,i),c(3,i)
163      end do
164
165      end
166
167      subroutine mm_open_file(filename,fn)
168      implicit none
169#include "util.fh"
170#include "errquit.fh"
171#include "inp.fh"
172#include "stdio.fh"
173      character*(*) filename
174      integer       fn
175c
176      character*180 buffer
177      character*180 message
178      character*30 pname,atag
179c
180      logical util_io_unit
181      external util_io_unit
182c
183      pname = "mm_open_file"
184c
185      if(.not.util_io_unit(80,90,fn))
186     +  call errquit(pname//"cannot get io unit",0,0)
187c     first try to open file in the run directory
188      buffer = filename
189      message = "opening file "//buffer
190      open(unit=fn,file=buffer,status='old',form="formatted",ERR=10)
191      goto 800
19210    continue
193c     now try perm directory
194      call util_file_name_resolve(buffer, .false.)
195      message = "opening file "//buffer
196      open(unit=fn,file=buffer,status='old',form="formatted",ERR=911)
197800   continue
198c     write(luout,*) "Successfully "//trim(message)
199c     write(luout,*)
200      return
201911   call errquit("error "//trim(message),0,
202     >        -1)
203      end
204
205      subroutine mm_add_energy(rtdb,e)
206      implicit none
207#include "util.fh"
208#include "errquit.fh"
209#include "inp.fh"
210#include "stdio.fh"
211#include "mafdecls.fh"
212#include "msgids.fh"
213#include "global.fh"
214#include "rtdb.fh"
215#include "mm_bq_data.fh"
216#include "mm_geom_data.fh"
217#include "bq.fh"
218
219      integer rtdb
220      double precision e
221
222      double precision bq_el_energy
223      double precision bq_nuc_energy0
224      double precision eqm0
225      double precision eqmtot
226      double precision emm
227      double precision ebq0
228c     dummy variables to pass as an argument in mm_cons_reaction
229      double precision g(3,nact)
230c
231      logical rtdb_mode
232      integer master
233      character*30 pname
234c
235      pname = "mm_add_energy"
236c     write(*,*) pname
237      master=0
238
239      if(.not.bq_nuc_energy(rtdb,ebq0))
240     >  call errquit('failed bq energy',0,RTDB_ERR)
241
242      call mm_geom_restore(rtdb)
243      call mm_cons_reaction(rtdb,.false.,e,g)
244
245      eqmtot = e
246      bq_nuc_energy0 = ebq0
247
248c     bq elec interaction energy
249      if (.not. rtdb_get(rtdb,'dft:bq_energy',mt_dbl,1,bq_el_energy))
250     $     bq_el_energy = 0.0d0
251
252      eqm0 = eqmtot - bq_el_energy - bq_nuc_energy0
253
254      rtdb_mode = rtdb_parallel(.false.)
255
256      if(ga_nodeid().eq.master) then
257        call mm_links_bq_add_energy(rtdb,e)
258c       call mm_links_bq_scaled_add_energy(rtdb,e)
259        call mm_bonded_add_energy(rtdb,e)
260        call mm_vdw_add_energy(rtdb,e)
261      end if
262
263      if(ga_nodeid().eq.master) then
264
265        emm = e - eqmtot
266
267        if(.not. rtdb_put(rtdb,'qmmm:qm_en',mt_dbl,1,eqmtot))
268     >  call errquit('qmmm: failed put qm energy',0,RTDB_ERR)
269
270        if(.not. rtdb_put(rtdb,'qmmm:qm_int_en',mt_dbl,1,eqm0))
271     >  call errquit('qmmm: failed put qm internal energy',0,RTDB_ERR)
272
273        if(.not. rtdb_put(rtdb,'qmmm:bq_nuc',mt_dbl,1,bq_nuc_energy0))
274     >  call errquit('qmmm: failed put bq_nuc energy',0,RTDB_ERR)
275
276        if(.not. rtdb_put(rtdb,'qmmm:bq_el',mt_dbl,1,bq_el_energy))
277     >  call errquit('qmmm: failed put bq_el energy',0,RTDB_ERR)
278
279        if(.not. rtdb_put(rtdb,'qmmm:mm_en',mt_dbl,1,emm))
280     >  call errquit('qmmm: failed put emm',0,RTDB_ERR)
281
282        if(.not. rtdb_put(rtdb,'qmmm:tot_en',mt_dbl,1,e))
283     >  call errquit('qmmm: failed put total energy',0,RTDB_ERR)
284
285        call mm_print_energy(rtdb)
286      end if
287
288      call ga_sync()
289      call ga_brdcst(msg_smd,e,
290     >     ma_sizeof(mt_dbl,1,mt_byte),master)
291      call ga_sync()
292      rtdb_mode = rtdb_parallel(rtdb_mode)
293
294      end
295
296      subroutine mm_add_egrad(rtdb,e,n,g)
297      implicit none
298#include "util.fh"
299#include "errquit.fh"
300#include "inp.fh"
301#include "stdio.fh"
302#include "mafdecls.fh"
303#include "msgids.fh"
304#include "global.fh"
305#include "rtdb.fh"
306#include "mm_coords_data.fh"
307#include "mm_link_data.fh"
308#include "mm_bq_data.fh"
309#include "bq.fh"
310
311      integer rtdb
312      double precision e
313
314      double precision bq_el_energy
315      double precision bq_nuc_energy0
316      double precision eqm0
317      double precision eqmtot
318      double precision emm
319      double precision ebq0
320
321      integer n
322      double precision g(3,n)
323c
324      integer i,j
325      character*30 pname
326      integer master
327      logical rtdb_mode
328      integer l_act, k_act
329      integer nactive
330      character*30 operation
331      logical ignore
332
333      pname = "mm_add_egrad"
334c     write(*,*) pname
335
336      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
337
338      master=0
339
340      if(ga_nodeid().eq.master) then
341        call mm_links_adjust_forces(n,int_mb(i_iqml),g)
342      end if
343      call ga_sync()
344
345      call mm_link_ebq_add_grad(rtdb,n,g)
346
347      if(.not.bq_nuc_energy(rtdb,ebq0))
348     >  call errquit('failed bq energy',0,RTDB_ERR)
349
350      bq_nuc_energy0 = ebq0
351
352      call mm_geom_restore(rtdb)
353      call mm_cons_reaction(rtdb,.true.,e,g)
354
355      eqmtot = e
356
357c     bq elec interaction energy
358      if (.not. rtdb_get(rtdb,'dft:bq_energy',mt_dbl,1,bq_el_energy))
359     $     bq_el_energy = 0.0d0
360
361      eqm0 = eqmtot - bq_el_energy - bq_nuc_energy0
362
363      rtdb_mode = rtdb_parallel(.false.)
364
365      if(ga_nodeid().eq.master) then
366C     lookup table and list of active atoms
367      if (.not. ma_push_get(MT_LOG,n,'active atoms',l_act,k_act))
368     $     call errquit('grad: could not allocate l_act',n, MA_ERR)
369        call mm_links_bq_add_egrad(rtdb,e,n,g)
370c       call mm_links_bq_scaled_add_egrad(rtdb,e,n,g)
371        call mm_bonded_add_grad(rtdb,e,n,g)
372        call mm_vdw_add_egrad(rtdb,e,n,g)
373
374        call grad_active_atoms(rtdb, n, log_mb(k_act), nactive)
375
376c       TP: callin zero_forces() for fixed atom
377c       defined in constraint module
378        call zero_forces(g,log_mb(k_act),n)
379
380        call mm_print_grad_tot(rtdb,n,g)
381
382      if (.not.ma_pop_stack(l_act))
383     $     call errquit('grad:ma free act',1, MA_ERR)
384      end if
385
386
387      if(ga_nodeid().eq.master) then
388
389        emm = e - eqmtot
390
391        if(.not. rtdb_put(rtdb,'qmmm:qm_en',mt_dbl,1,eqmtot))
392     >  call errquit('qmmm: failed put qm energy',0,RTDB_ERR)
393
394        if(.not. rtdb_put(rtdb,'qmmm:qm_int_en',mt_dbl,1,eqm0))
395     >  call errquit('qmmm: failed put qm internal energy',0,RTDB_ERR)
396
397        if(.not. rtdb_put(rtdb,'qmmm:bq_nuc',mt_dbl,1,bq_nuc_energy0))
398     >  call errquit('qmmm: failed put bq_nuc energy',0,RTDB_ERR)
399
400        if(.not. rtdb_put(rtdb,'qmmm:bq_el',mt_dbl,1,bq_el_energy))
401     >  call errquit('qmmm: failed put bq_el energy',0,RTDB_ERR)
402
403        if(.not. rtdb_put(rtdb,'qmmm:mm_en',mt_dbl,1,emm))
404     >  call errquit('qmmm: failed put emm',0,RTDB_ERR)
405
406        if(.not. rtdb_put(rtdb,'qmmm:tot_en',mt_dbl,1,e))
407     >  call errquit('qmmm: failed put total energy',0,RTDB_ERR)
408
409        call mm_print_energy(rtdb)
410        call mm_create_xyz_file(rtdb)
411C       TP: create intermediate crdparms file for each bead for next iter
412        if(operation.eq."neb") call mm_create_crdparm(rtdb)
413      end if
414
415      call ga_sync()
416      call ga_brdcst(msg_smd,e,
417     >     ma_sizeof(mt_dbl,1,mt_byte),master)
418      call ga_brdcst(msg_smd,g,
419     >     3*n*ma_sizeof(mt_dbl,1,mt_byte),master)
420      call ga_sync()
421
422      rtdb_mode = rtdb_parallel(rtdb_mode)
423
424c     call mm_create_xyz_file(rtdb)
425
426      return
427911   call errquit("error "//trim(pname),0,-1)
428      end
429
430      subroutine mm_bonded_add_energy(rtdb,e)
431      implicit none
432#include "util.fh"
433#include "errquit.fh"
434#include "inp.fh"
435#include "stdio.fh"
436#include "mafdecls.fh"
437#include "msgids.fh"
438#include "global.fh"
439#include "rtdb.fh"
440      integer rtdb
441      double precision e
442      character*30 pname
443
444      pname = "mm_bonded_add_energy"
445c     write(*,*) pname
446
447      call mm_bond_add_energy(rtdb,e)
448      call mm_angl_add_energy(rtdb,e)
449      call mm_dihe_add_energy(rtdb,e)
450
451      end
452
453
454      subroutine mm_bonded_add_grad(rtdb,e,n,g)
455      implicit none
456#include "util.fh"
457#include "errquit.fh"
458#include "inp.fh"
459#include "stdio.fh"
460#include "mafdecls.fh"
461#include "msgids.fh"
462#include "global.fh"
463#include "rtdb.fh"
464      integer rtdb
465      double precision e
466      integer n
467      double precision g(n)
468      character*30 pname
469
470      pname = "mm_bonded_add_grad"
471c     write(*,*) pname
472
473      call mm_bond_add_egrad(rtdb,e,n,g)
474      call mm_angl_add_egrad(rtdb,e,n,g)
475      call mm_dihe_add_egrad(rtdb,e,n,g)
476
477      end
478
479      subroutine mm_print_grad_tot(rtdb,n,g)
480      implicit none
481#include "util.fh"
482#include "errquit.fh"
483#include "inp.fh"
484#include "stdio.fh"
485#include "mafdecls.fh"
486#include "geom.fh"
487#include "rtdb.fh"
488#include "mm_geom_data.fh"
489#include "mm_link_data.fh"
490#include "mm_coords_data.fh"
491      integer rtdb
492      integer n
493      double precision g(3,n)
494
495      logical status
496      logical ignore
497      integer geom
498      integer i_c,h_c
499      integer i_tag,h_tag
500      integer i, j, ncent
501      character*30 pname
502      character*30 message
503      character*16 tag
504      double precision scale
505      logical geom_cart_get1
506      external geom_cart_get1
507
508      pname="mm_print_grad_tot"
509c     write(*,*) pname
510
511      if(.not.ma_push_get(mt_dbl,3*nact,'c',h_c,i_c))
512     & call errquit('mm: Failed to allocate memory for c',
513     & 3*nact, MA_ERR)
514
515      call dcopy(nact*3,dbl_mb(i_rqml),1,dbl_mb(i_c),1)
516      call util_convert_units("angstrom","au",scale)
517      call dscal(3*nact, scale,dbl_mb(i_c),1)
518
519      write(6,1000) "QM + MM",
520     $     'x','y','z','x','y','z'
521      do i = 1, nact
522         call mm_coords_tag_set(byte_mb(i_tqml+16*(i-1)),tag)
523         write(6,2000) i, tag,
524     $        (dbl_mb(i_c+3*(i-1)+j),j=0,2),
525     $         g(1,i),g(2,i),g(3,i)
526      enddo
527      write(6,*)
528
529 1000 format(/,/,25X,A,' ENERGY GRADIENTS',/,/,4X,'atom',15X,
530     $     'coordinates',
531     $     24X,'gradient',/,6X,2(1X,(3(10X,A1))))
532 2000 format(1X,I3,1X,A4,2(1X,3(1X,F10.6)))
533      write(6,*)
534      call util_flush(6)
535
536      message = "memory deallocation"
537      if(.not.ma_pop_stack(h_c)) goto 911
538
539      return
540
541911   call errquit("error "//trim(message),0,-1)
542
543      end
544
545      subroutine mm_print_energy(rtdb)
546      implicit none
547#include "util.fh"
548#include "mafdecls.fh"
549#include "errquit.fh"
550#include "rtdb.fh"
551#include "msgids.fh"
552#include "global.fh"
553#include "stdio.fh"
554
555      integer rtdb
556
557      character*32 pname
558
559      double precision bq_el_energy
560      double precision bq_nuc_energy
561      double precision bq_energy
562      double precision eqm0
563      double precision eqmtot
564      double precision emm
565      double precision etot
566
567      double precision scale_energy
568
569      pname = "mm_print_energy"
570c     write(*,*) pname
571
572      if(ga_nodeid().ne.0) return
573
574      if(.not. rtdb_get(rtdb,'qmmm:qm_en',mt_dbl,1,eqmtot))
575     >  call errquit('qmmm: failed get qm energy',0,RTDB_ERR)
576
577      if(.not. rtdb_get(rtdb,'qmmm:qm_int_en',mt_dbl,1,eqm0))
578     >  call errquit('qmmm: failed get qm internal energy',0,RTDB_ERR)
579
580      if(.not. rtdb_get(rtdb,'qmmm:bq_nuc',mt_dbl,1,bq_nuc_energy))
581     >  call errquit('qmmm: failed get bq_nuc energy',0,RTDB_ERR)
582
583      if(.not. rtdb_get(rtdb,'qmmm:bq_el',mt_dbl,1,bq_el_energy))
584     >  call errquit('qmmm: failed get bq_el energy',0,RTDB_ERR)
585
586      if(.not. rtdb_get(rtdb,'qmmm:mm_en',mt_dbl,1,emm))
587     >  call errquit('qmmm: failed get emm',0,RTDB_ERR)
588
589      if(.not. rtdb_get(rtdb,'qmmm:tot_en',mt_dbl,1,etot))
590     >  call errquit('qmmm: failed get total energy',0,RTDB_ERR)
591
592      bq_energy = bq_el_energy + bq_nuc_energy
593
594      call util_convert_units("au","kcal",scale_energy)
595
596      write(*,19)
597      write(*,21) "QM/MM Energy"
598      write(*,19)
599      write(*,21) "quantum energy", eqmtot, eqmtot*scale_energy
600c     write(*,21) "quantum energy adjusted", qm_energy-eatoms,
601c    >            (qm_energy-eatoms)*cau2kj
602      if(bq_energy.ne.0.0d0) then
603        write(*,21) "quantum energy internal", eqm0,
604     >               eqm0*scale_energy
605        write(*,21) "Bq-nuclear energy", bq_nuc_energy,
606     >               bq_nuc_energy*scale_energy
607        write(*,21) "Bq-electron energy", bq_el_energy,
608     >               bq_el_energy*scale_energy
609      end if
610      write(*,21) "classical energy", emm,
611     >            emm*scale_energy
612      write(*,21) "total qmmm energy", etot, etot*scale_energy
613      write(*,19)
614
615      write(*,*)
61619    FORMAT("@",72("-"))
61720    FORMAT("@",1X,A,T34,F18.9)
61821    FORMAT("@",1X,A,T34,F18.9," (",E12.6,2X, "kcal/mol)")
619
620      end
621
622      subroutine mm_create_xyz_file(rtdb)
623      implicit none
624#include "util.fh"
625#include "mafdecls.fh"
626#include "errquit.fh"
627#include "rtdb.fh"
628#include "msgids.fh"
629#include "global.fh"
630#include "stdio.fh"
631      integer rtdb
632
633      character*50 filename
634      character*255 full_filename
635      character*30 operation
636      logical ignore
637      character*32 pname
638      character*20 tag
639      integer ibead
640      logical wrt_neb_xyz
641      logical wrt_opt_xyz
642
643      pname = "mm_create_xyz_file"
644c     write(*,*) pname
645
646      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
647
648      wrt_neb_xyz = .false.
649      wrt_opt_xyz = .false.
650
651      if(operation.eq."neb") wrt_neb_xyz = .true.
652      if(operation.eq."optimize") wrt_opt_xyz = .true.
653      if(operation.eq."saddle") wrt_opt_xyz = .true.
654
655      if(wrt_neb_xyz) then
656        if(.not.rtdb_get(rtdb,'neb:ibead',mt_int,1,ibead))
657     >    call errquit('neb:ibead get',0,-1)
658        tag = ' '
659        write(tag,10) ibead
660        call util_file_prefix(tag,filename)
661        call util_file_name_noprefix(filename,.false.,
662     >                               .false.,
663     >                               full_filename)
664
665        open(unit=19,file=full_filename,form='formatted')
666        call mm_neb_create_xyz_file(19,ibead)
667        close(19)
668
669      else if(wrt_opt_xyz) then
670        call mm_geom_create_xyz_file(rtdb)
671      else
672        return
673      end if
674
67510    format('neb_',i3.3,'.xyz')
676
677      end
678
679      subroutine mm_create_crdparm(rtdb)
680      implicit none
681#include "util.fh"
682#include "mafdecls.fh"
683#include "errquit.fh"
684#include "rtdb.fh"
685#include "msgids.fh"
686#include "global.fh"
687#include "stdio.fh"
688      integer rtdb
689
690      character*50 filename
691      character*255 full_filename
692      character*30 operation
693      logical ignore
694      character*32 pname
695      character*20 tag
696      integer ibead
697
698      pname = "mm_create_crdparm"
699c     write(*,*) pname
700
701      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
702
703      if(operation.eq."neb") then
704        if(.not.rtdb_get(rtdb,'neb:ibead',mt_int,1,ibead))
705     >    call errquit('neb:ibead get',0,-1)
706        filename = ' '
707        write(filename,10) ibead
708
709
710      else if(operation.eq."optimize".or.operation.eq."saddle") then
711        if(.not.rtdb_cget(rtdb,"mm:crdparms:load",1,
712     &                   filename))
713     &  call errquit("error "//"mm:crdparmfile",0,-1)
714
715      else
716        return
717      end if
718
719      open(unit=19,file=filename,form='formatted')
720      call mm_write_qmcoords(19)
721      call mm_write_mmcoords(19)
722      call mm_write_bond(19)
723      call mm_write_angl(19)
724      call mm_write_dihe(19)
725      call mm_write_vdw(19)
726      call mm_write_vdw14(19)
727      close(19)
728
72910    format('nwchem_',i3.3,'.crdparms')
730
731      end
732
733      subroutine mm_task_gradient(rtdb,theory,e,g,status)
734#include "util.fh"
735#include "mafdecls.fh"
736#include "errquit.fh"
737#include "rtdb.fh"
738#include "msgids.fh"
739#include "global.fh"
740#include "stdio.fh"
741
742c     input variables
743      integer rtdb
744      character*32 theory
745      double precision e
746      double precision g(3,*)
747      logical status
748
749c     local variables
750      logical ignore
751      character*30 operation
752
753      logical task_gradient_doit
754      external task_gradient_doit
755
756      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
757      if(operation.eq."neb")
758     >   call mm_neb_update_ibead(rtdb)
759      call mm_geom_push_active(rtdb)
760      call mm_geom_create_full(rtdb)
761      call mm_map_fixed_atoms(rtdb)
762      call mm_link_update_bq_coords(rtdb)
763      status = task_gradient_doit(rtdb,theory,e,g)
764      call mm_restore_fixed_atoms(rtdb)
765
766      end
767
768      subroutine mm_task_energy(rtdb,theory,e,status)
769#include "util.fh"
770#include "mafdecls.fh"
771#include "errquit.fh"
772#include "rtdb.fh"
773#include "msgids.fh"
774#include "global.fh"
775#include "stdio.fh"
776
777c     input variables
778      integer rtdb
779      character*32 theory
780      double precision e
781      logical status
782
783c     local variables
784      logical ignore
785      character*30 operation
786
787      logical task_energy_doit
788      external task_energy_doit
789
790      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
791      if(operation.eq."neb")
792     >   call mm_neb_update_ibead(rtdb)
793      call mm_geom_push_active(rtdb)
794      call mm_geom_create_full(rtdb)
795      call mm_link_update_bq_coords(rtdb)
796      status = task_energy_doit(rtdb,theory,e)
797
798      end
799