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