1 2 subroutine mm_geom_init(rtdb,geomname) 3 implicit none 4#include "mafdecls.fh" 5#include "errquit.fh" 6#include "geom.fh" 7#include "rtdb.fh" 8#include "util.fh" 9#include "global.fh" 10#include "inp.fh" 11#include "mm_coords_data.fh" 12#include "mm_bond_coords_data.fh" 13#include "mm_link_data.fh" 14#include "mm_geom_data.fh" 15 16 integer rtdb 17 character*(*) geomname 18 19 character*255 aname 20 character*32 pname 21 character*30 operation 22 logical orestrt 23 integer geom 24 integer i, j 25 logical ignore 26 integer i_c,h_c 27 integer i_tag,h_tag 28 integer i_m,h_m 29 integer i_atn,h_atn 30 integer i_dbq,h_dbq 31 double precision scale 32 logical geom_tag_to_atn 33 logical geom_tag_to_charge 34 external geom_tag_to_atn 35 external geom_tag_to_charge 36 37 38 pname = "mm_geom_init" 39c write(*,*) pname 40 41 ignore = rtdb_cget(rtdb,'task:operation',1,operation) 42 43c deallocate all previous allocated arrays just in case 44 call mm_geom_end() 45 46 nact = nqm 47 if(qmlink) nact = nqm + nlink 48 49 nfg = nqm + nlink 50 51 aux_geom = .false. 52c TP: refer to qmmm.F for operation = "hessian" and 53c operation = "freq". aux_geom = .true. 54 if(operation.eq."optimize") aux_geom =.true. 55 if(operation.eq."neb") aux_geom =.true. 56 if(operation.eq."freq") aux_geom =.true. 57 if(operation.eq."hessian") aux_geom =.true. 58 if(operation.eq."saddle") aux_geom =.true. 59 60c initialize indexing for a full qm geometry 61 if(.not.ma_alloc_get(mt_int,nact,'mm act atom ind', 62 & h_iact,i_iact)) 63 & call errquit(pname//'Failed to allocate heap',nact, 64 & MA_ERR) 65 66 if(.not.ma_alloc_get(mt_int,nfg,'mm fullg ind',h_ifg,i_ifg)) 67 & call errquit(pname//'Failed to allocate heap',nfg, 68 & MA_ERR) 69 70 call icopy(nqm,int_mb(i_iqm),1,int_mb(i_ifg),1) 71 call icopy(nlink,int_mb(i_lb+nlink),1,int_mb(i_ifg+nqm),1) 72 73 do i=1,nact 74 int_mb(i_iact+i-1) = int_mb(i_iqml+i-1) 75 end do 76 77 if(operation.eq."neb") then 78 if(.not.rtdb_get(rtdb,"mm:neb:restart",mt_log,1,orestrt)) 79 > orestrt = .false. 80 if(orestrt) return 81 end if 82 83 84 if(.not.ma_push_get(mt_dbl,3*nact,'c',h_c,i_c)) 85 & call errquit('mm: Failed to allocate memory for c', 86 & 3*nact, MA_ERR) 87 if(.not.ma_push_get(mt_dbl,nact,'q',h_dbq,i_dbq)) 88 & call errquit('mm: Failed to allocate memory for q',nact, 89 & MA_ERR) 90 if(.not.ma_push_get(mt_dbl,nact,'m',h_m,i_m)) 91 & call errquit('mm: Failed to allocate memory for m',nact, 92 & MA_ERR) 93 if(.not.ma_push_get(mt_int,nact,'inum',h_atn,i_atn)) 94 & call errquit('mm: Failed to allocate memory for atn',nact, 95 & MA_ERR) 96 if(.not.ma_push_get(mt_byte,16*nact,'t',h_tag,i_tag)) 97 & call errquit('mm: Failed to allocate memory for t',nact, 98 & MA_ERR) 99 100c assign coordinates for active atoms 101 do i=1,nact 102 dbl_mb(i_c+(i-1)*3) = dbl_mb(i_rqml+3*(i-1)) 103 dbl_mb(i_c+(i-1)*3+1) = dbl_mb(i_rqml+3*(i-1)+1) 104 dbl_mb(i_c+(i-1)*3+2) = dbl_mb(i_rqml+3*(i-1)+2) 105 call mm_coords_tag_set(byte_mb(i_tqml+16*(i-1)), 106 & byte_mb(i_tag+16*(i-1))) 107 end do 108 109 call util_convert_units("angstrom","au",scale) 110 call dscal(3*nact, scale,dbl_mb(i_c),1) 111 112 if(.not.geom_tag_to_charge(nact,byte_mb(i_tag),dbl_mb(i_dbq))) 113 & call errquit('mm: Failed to get charge from tag',nact,MA_ERR) 114 115 if(.not.geom_tag_to_atn(nact,byte_mb(i_tag),int_mb(i_atn))) 116 & call errquit('mm: Failed to get atn from tag',nact,MA_ERR) 117 118c if nlink > 0 assign empirical charge to mmlink 119 if(qmlink) then 120 do i=1,nlink 121 dbl_mb(i_dbq+nqm+i-1) = dbl_mb(i_lnkchg+i-1) 122 end do 123 end if 124 125 call mm_atn_get_mass(nact,int_mb(i_atn),dbl_mb(i_m)) 126 127c ignore = rtdb_delete(rtdb,"geometry") 128 ignore = rtdb_delete(rtdb,geomname(1:inp_strlen(geomname))) 129 130c if(.not.geom_create(geom,"geometry")) 131 if(.not.geom_create(geom,geomname(1:inp_strlen(geomname)))) 132 & call errquit('mm: Failed to create geometry',0, GEOM_ERR) 133 134 if(.not.geom_cart_set(geom,nact,byte_mb(i_tag), 135 & dbl_mb(i_c), 136 & dbl_mb(i_dbq))) 137 & call errquit('mm: Failed to initialize geometry',0, GEOM_ERR) 138 139 if(.not.geom_masses_set(geom,nact,dbl_mb(i_m))) 140 & call errquit('mm: Failed to initialize masses',0, GEOM_ERR) 141 call geom_compute_values(geom) 142 143 if(.not.geom_rtdb_store(rtdb,geom, 144 & geomname(1:inp_strlen(geomname)))) 145 & call errquit('mm: Failed to store geom to rtdb',0, RTDB_ERR) 146 147 if(.not.geom_destroy(geom)) 148 & call errquit('mm: Failed to destroy geometry',0, GEOM_ERR) 149 150 if(.not.ma_pop_stack(h_tag)) 151 & call errquit('mm: Failed to deallocate stack t_all',nact, 152 & MA_ERR) 153 if(.not.ma_pop_stack(h_atn)) 154 & call errquit('mm: Failed to deallocate stack atn_all',nact, 155 & MA_ERR) 156 if(.not.ma_pop_stack(h_m)) 157 & call errquit('mm: Failed to deallocate stack m_all',nact, 158 & MA_ERR) 159 if(.not.ma_pop_stack(h_dbq)) 160 & call errquit('mm: Failed to deallocate stack q_all',nact, 161 & MA_ERR) 162 if(.not.ma_pop_stack(h_c)) 163 & call errquit('mm: Failed to deallocate stack c_all',nact, 164 & MA_ERR) 165 166 end 167 168 subroutine mm_geom_end() 169 implicit none 170#include "mafdecls.fh" 171#include "errquit.fh" 172#include "geom.fh" 173#include "rtdb.fh" 174#include "global.fh" 175#include "inp.fh" 176#include "mm_geom_data.fh" 177 178 character*30 pname 179 pname = "mm_geom_end" 180 181 if(nfg.gt.0) then 182 if (.not.ma_free_heap(h_ifg)) goto 911 183 if (.not.ma_free_heap(h_iact)) goto 911 184 nfg = 0 185 nact = 0 186 aux_geom = .false. 187 end if 188 189 return 190911 call errquit("error "//trim(pname),0,-1) 191 192 end 193 194 subroutine mm_geom_create_full(rtdb) 195 implicit none 196#include "mafdecls.fh" 197#include "errquit.fh" 198#include "geom.fh" 199#include "rtdb.fh" 200#include "global.fh" 201#include "inp.fh" 202#include "mm_geom_data.fh" 203 204 character*32 pname 205 integer rtdb 206 integer i_c,h_c 207 integer i_t,h_t 208 integer i_m,h_m 209 integer i_atn,h_atn 210 integer i_dbq,h_dbq 211 212 pname = "geom_create_full" 213c write(*,*) 'in ', pname 214 215 if(.not.ma_push_get(mt_dbl,3*nfg,'c',h_c,i_c)) 216 & call errquit('mm: Failed to allocate memory for c', 217 & 3*nfg, MA_ERR) 218 if(.not.ma_push_get(mt_dbl,nfg,'q',h_dbq,i_dbq)) 219 & call errquit('mm: Failed to allocate memory for q',nfg, 220 & MA_ERR) 221 if(.not.ma_push_get(mt_dbl,nfg,'m',h_m,i_m)) 222 & call errquit('mm: Failed to allocate memory for m',nfg, 223 & MA_ERR) 224 if(.not.ma_push_get(mt_int,nfg,'inum',h_atn,i_atn)) 225 & call errquit('mm: Failed to allocate memory for atn',nfg, 226 & MA_ERR) 227 if(.not.ma_push_get(mt_byte,16*nfg,'t',h_t,i_t)) 228 & call errquit('mm: Failed to allocate memory for t',nfg, 229 & MA_ERR) 230 231 call mm_get_geom(rtdb,nfg,int_mb(i_ifg),dbl_mb(i_c), 232 & dbl_mb(i_dbq),dbl_mb(i_m), 233 & int_mb(i_atn),byte_mb(i_t)) 234 235 if(.not.rtdb_cget(rtdb,'geometry',1,oldgeom)) 236 & oldgeom = ' ' 237 238 if(.not.rtdb_cput(rtdb,'geometry',1,'full_geom')) 239 & call errquit(pname//' storing geom name to rtdb',0, RTDB_ERR) 240 241c release temporary memory 242c ------------------------ 243 if(.not.ma_pop_stack(h_t)) 244 & call errquit('mm: Failed to deallocate stack t_all',nfg, 245 & MA_ERR) 246 if(.not.ma_pop_stack(h_atn)) 247 & call errquit('mm: Failed to deallocate stack atn_all',nfg, 248 & MA_ERR) 249 if(.not.ma_pop_stack(h_m)) 250 & call errquit('mm: Failed to deallocate stack m_all',nfg, 251 & MA_ERR) 252 if(.not.ma_pop_stack(h_dbq)) 253 & call errquit('mm: Failed to deallocate stack q_all',nfg, 254 & MA_ERR) 255 if(.not.ma_pop_stack(h_c)) 256 & call errquit('mm: Failed to deallocate stack c_all',nfg, 257 & MA_ERR) 258 259 260 end 261 262 subroutine mm_get_geom(rtdb,nt,ifg,c,dbq,m,atn,t) 263 implicit none 264#include "mafdecls.fh" 265#include "errquit.fh" 266#include "geom.fh" 267#include "rtdb.fh" 268#include "util.fh" 269#include "mm_bond_coords_data.fh" 270#include "mm_coords_data.fh" 271 272 integer rtdb 273 integer geom 274 integer nt 275 integer ifg(nt) 276 double precision c(3,nt) 277 double precision dbq(nt) 278 double precision m(nt) 279 integer atn(nt) 280 character*16 t(nt) 281 282 integer i 283 integer nat 284 double precision scale 285 logical status 286 logical ignore 287 character*30 message 288 character*30 pname 289 290 pname = "mm_get_geom" 291 292c if (.not. geom_create(geom, 'full_geom')) 293 if (.not. geom_create(geom, 'geometry')) 294 & call errquit('cons_create_geom',0, GEOM_ERR) 295 296c if (.not. geom_rtdb_load(rtdb, geom, 'full_geom')) 297 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 298 & call errquit('cons_load_geom',0, RTDB_ERR) 299 300 status=geom_ncent(geom,nat) 301 302 if(.not.status) 303 & call errquit('cons_init: geom_create?',70, GEOM_ERR) 304 305 if(.not.geom_cart_get2(geom,nat,t,c,dbq,atn)) 306 & goto 911 307 308 call util_convert_units("au","angstrom",scale) 309 call dscal(3*nat, scale,c,1) 310 311 if(.not.geom_destroy(geom)) 312 & goto 911 313 314 ignore = geom_rtdb_delete(rtdb,"full_geom") 315 316 if(.not.geom_create(geom,"full_geom")) 317 & call errquit('mm: Failed to create geometry',0, GEOM_ERR) 318 319 call mm_links_adjust(nt,ifg,atn,t,c,dbq) 320 321 322 call util_convert_units("angstrom","au",scale) 323 call dscal(3*nt, scale,c,1) 324 325 if(.not.geom_cart_set(geom,nt,t,c,dbq)) 326 & call errquit('mm: Failed to initialize geometry',0, GEOM_ERR) 327 328 call mm_atn_get_mass(nt,atn,m) 329 330 if(.not.geom_masses_set(geom,nt,m)) 331 & call errquit('mm: Failed to initialize masses',0, GEOM_ERR) 332 call geom_compute_values(geom) 333 334 if(.not.geom_rtdb_store(rtdb,geom,"full_geom")) 335 & call errquit('mm: Failed to store geom to rtdb',0, RTDB_ERR) 336 337 if(.not.geom_destroy(geom)) 338 & call errquit('mm: Failed to destroy geometry',0, GEOM_ERR) 339 340 return 341911 call errquit("error "//trim(message),0,-1) 342 end 343 344 subroutine mm_atn_get_mass(n,atn,m) 345 implicit none 346#include "mafdecls.fh" 347#include "errquit.fh" 348#include "geom.fh" 349 integer n 350 integer atn(n) 351 double precision m(n) 352 353c local variables 354 integer i 355 character*32 pname 356 357 pname = "mm_atn_get_mass" 358 359 do i=1,n 360 if(.not.geom_atn_to_default_mass(atn(i),m(i))) 361 & call errquit(pname,0, GEOM_ERR) 362 363 end do 364 365 end 366 367 subroutine mm_geom_push_active(rtdb) 368 implicit none 369#include "mafdecls.fh" 370#include "errquit.fh" 371#include "geom.fh" 372#include "rtdb.fh" 373#include "global.fh" 374#include "inp.fh" 375#include "mm_geom_data.fh" 376#include "mm_bond_coords_data.fh" 377#include "mm_coords_data.fh" 378 379 integer rtdb 380 381 integer geom 382 integer ncent 383 integer i_ctmp,h_ctmp 384 integer i, j 385 integer nt 386 double precision scale 387 character*32 pname 388 389 pname = "mm_geom_push_active" 390c write(*,*) 'in ', pname 391 392 nt = nact 393 394c -------------------------------------- 395c get active coordinates out of geometry 396c -------------------------------------- 397 if(.not.geom_create(geom,'geometry')) 398 & call errquit('mm: Failed to create geometry',0, GEOM_ERR) 399 400 if(.not.geom_rtdb_load(rtdb,geom,"geometry")) 401 & call errquit('mm: Failed to create geometry',0, GEOM_ERR) 402 403 if(.not. geom_ncent(geom, ncent) ) 404 & call errquit("mm:geom_ncent",0,0) 405 406 if(.not.ma_push_get(mt_dbl,3*ncent,'ctmp',h_ctmp,i_ctmp)) 407 & call errquit( pname//'Failed to allocate memory for ctmp', 408 & 3*ncent, MA_ERR) 409 410 if(.not. geom_cart_coords_get(geom,dbl_mb(i_ctmp))) 411 & call errquit("mm:geom_cart_coords_get",0,0) 412 413 414 call util_convert_units("au","angstrom",scale) 415 call dscal(3*ncent, scale,dbl_mb(i_ctmp),1) 416 call mm_set_coord(nt,int_mb(i_iact),dbl_mb(i_ctmp)) 417 call mm_vdw_qmcoords_load(rtdb) 418 call mm_vdw14_qmcoords_load(rtdb) 419 call mm_bond_qmcoords_load(rtdb) 420 421 if(.not.ma_pop_stack(h_ctmp)) 422 & call errquit('mm: 423 & Failed to deallocate stack c_tmp',ncent, 424 & MA_ERR) 425 426 if(.not.geom_destroy(geom)) 427 & call errquit('mm: Failed to destroy geometry',0, GEOM_ERR) 428 429 430 end 431 432 subroutine mm_set_coord(nt,ai,c) 433 implicit none 434#include "mafdecls.fh" 435#include "errquit.fh" 436#include "geom.fh" 437#include "rtdb.fh" 438#include "global.fh" 439#include "inp.fh" 440#include "mm_bond_coords_data.fh" 441#include "mm_link_data.fh" 442#include "mm_geom_data.fh" 443#include "mm_coords_data.fh" 444 445 integer rtdb 446 integer nt 447 integer ai(nt) 448 double precision c(3,nt) 449 450 integer i,j,k 451 integer indx, iqm, imm 452 character*32 pname 453 454 pname = "mm_set_coord" 455c write(*,*) 'in ', pname 456 457 do i=1,nt 458 do j=1,nqml 459 indx = int_mb(i_iqml+j-1) 460 if(ai(i).eq.indx) then 461 do k = 1,3 462 dbl_mb(i_rqml+3*(j-1)+k-1) = c(k,i) 463 end do 464 end if 465 end do 466 end do 467 468c ------------------- 469c update i_rqm array 470c ------------------- 471 472 do i=1,nqm 473 iqm = int_mb(i_iqm+i-1) 474 do j=1,nqml 475 indx = int_mb(i_iqml+j-1) 476 if(iqm.eq.indx) then 477 do k = 1,3 478 dbl_mb(i_rqm+3*(i-1)+k-1) = dbl_mb(i_rqml+3*(j-1)+k-1) 479 end do 480 exit 481 end if 482 end do 483 end do 484 485c ------------------- 486c update i_rmm array 487c ------------------- 488 if(qmlink) then 489 do i=1,nlink 490 indx = int_mb(i_lb+nlink+i-1) 491 do j=1,nmm 492 imm = int_mb(i_imm+j-1) 493 if(imm.eq.indx) then 494 do k = 1,3 495 dbl_mb(i_rmm+3*(j-1)+k-1) = dbl_mb(i_rqml+ 496 > (nqm+i-1)*3+k-1) 497 end do 498 exit 499 end if 500 end do 501 end do 502 end if 503 504 end 505 506 subroutine mm_geom_restore(rtdb) 507 implicit none 508#include "mafdecls.fh" 509#include "errquit.fh" 510#include "geom.fh" 511#include "rtdb.fh" 512#include "mm_geom_data.fh" 513 integer rtdb 514 logical ignore 515 character*32 pname 516 517 pname = "mm_restore_geom" 518c write(*,*) pname 519 520 if(.not.aux_geom) return 521 522 ignore = rtdb_delete(rtdb,'geometry') 523 524 if(oldgeom.ne.' ') then 525 if(.not.rtdb_cput(rtdb,'geometry',1,oldgeom)) 526 & call errquit(pname//' storing geom to rtdb', 0, RTDB_ERR) 527 end if 528 529 end 530 531 subroutine mm_geom_create_xyz_file(rtdb) 532 implicit none 533#include "util.fh" 534#include "mafdecls.fh" 535#include "errquit.fh" 536#include "rtdb.fh" 537#include "msgids.fh" 538#include "global.fh" 539#include "stdio.fh" 540#include "geom.fh" 541 542 integer rtdb 543 544 integer geom 545 character*30 pname 546 character*30 geomname 547 character*30 operation 548 character*50 filename 549 character*255 full_filename 550 integer ncent 551 logical status 552 logical ignore 553 554 pname = "mm_geom_create_xyz_file" 555 geomname = "geometry" 556 557 558 if (ga_nodeid().eq.0) then 559 call util_file_prefix('opt.xyz',filename) 560 call util_file_name_noprefix(filename,.false., 561 > .false., 562 > full_filename) 563 564 if(.not.geom_create(geom,'geom_tmp')) 565 > call errquit('mm: Failed to create geometry',0, GEOM_ERR) 566 if(.not.geom_rtdb_load(rtdb,geom,geomname)) 567 > call errquit('mm: Failed to load geometry',0, GEOM_ERR) 568 status = geom_ncent(geom,ncent) 569 open(88,file=full_filename,form='formatted') 570 if (.not. geom_print_xyz(geom, 88)) 571 $ call errquit('mm:geom_print_xyz?',0, GEOM_ERR) 572 close(88) 573 574 if(.not.geom_destroy(geom)) 575 > call errquit('mm: Failed to destroy geometry',0, GEOM_ERR) 576 577 endif 578 579 end 580 581 subroutine mm_map_fixed_atoms(rtdb) 582 implicit none 583#include "rtdb.fh" 584#include "util.fh" 585#include "inp.fh" 586#include "mafdecls.fh" 587#include "errquit.fh" 588#include "geom.fh" 589#include "mm_link_data.fh" 590#include "mm_geom_data.fh" 591 592 integer rtdb 593 594 character*30 pname 595 integer i,j,k 596 integer nact0, nact1 597 integer i_cons,h_cons 598 integer i_act,h_act 599 integer ma_type 600 integer jj0,jj,ii,ii0 601 integer nal,h_al,i_al 602 integer h_am,i_am 603 logical aflag 604 605 pname = "mm_map_fixed_atoms" 606c write(*,*) pname 607 608c --------------- 609c get active list 610c --------------- 611 if (rtdb_ma_get(rtdb, 'geometry:actlist', ma_type, 612 $ nact0, h_cons)) then 613 if (.not. ma_get_index(h_cons, i_cons)) 614 $ call errquit(pname,h_cons, 615 & MA_ERR) 616 617 if (.not.rtdb_delete(rtdb, 'geometry:actlist')) 618 $ call errquit(pname,0, 619 & RTDB_ERR) 620 do i=1,nact0 621 end do 622 if (.not.rtdb_put(rtdb, 'qmmm:actlist', 623 > mt_int,nact0,int_mb(i_cons))) 624 $ call errquit(pname,0, 625 & RTDB_ERR) 626 627 else 628 nact0 = nact 629 if(.not.ma_alloc_get(mt_int, nact0, 'qmmm actlist', 630 & h_cons, i_cons) ) call errquit( 631 & 'qmmm_data_alloc: unable to allocate heap space', 632 & nact, MA_ERR) 633 do i=1,nact0 634 int_mb(i_cons+i-1) = i 635 end do 636 end if 637 638c ---------------------- 639c create new active list 640c ---------------------- 641 if(.not.ma_push_get(MT_INT, nfg, 'qmmm fixed atom list', 642 & h_act, i_act) ) call errquit( 643 & pname//' unable to allocate stack', 644 & nact, MA_ERR) 645 646c if(.not.ma_push_get(MT_LOG, ng, 'tmp qmmm act atom list', 647c & h_am, i_am) ) call errquit( 648c & pname//' unable to allocate stack', 649c & ng, MA_ERR) 650 651c TP: call following subroutine to map active atom. 652c For now, assume all (qm) atoms are active. 653 654c call mm_cons_get_map(ng,int_mb(i_ig), 655c > log_mb(i_am)) 656 657c -------------------------------------------- 658c find total number (nal) of link atoms 659c that are bonded to active qm atoms 660c and store their global index in i_al 661c Note that link atoms are always indexed 662c by the global index of the classical atom 663c -------------------------------------------- 664 if(nlink.ne.0) then 665 666 if(.not.ma_push_get(MT_INT, nlink, 'qmmm tmp link', 667 & h_al, i_al) ) call errquit( 668 & pname//' unable to allocate stack', 669 & nact, MA_ERR) 670 671 nal=0 672 do i=1,nact0 673c get global index next active atom in the auxilary geometry 674 ii0 = int_mb(i_cons+i-1) 675 ii = int_mb(i_iqml+ii0-1) 676 aflag = .true. 677 if(aflag) then 678 do j=1,nlink 679c get index of qm atom bonded to a link atom 680 jj = int_mb(i_lb+j-1) 681c if qm atom bonded to a link atom is active 682c store corresponding index of the link atom 683 if(aflag.and.(jj.eq.ii)) then 684 nal=nal+1 685 int_mb(i_al+nal-1)=int_mb(i_lb+nlink+j-1) 686 end if 687 end do 688 end if 689 end do 690 else 691 nal = 0 692 end if 693 694c ------------------------------------------------ 695 696c ------------------------------------------------ 697c construct active atom index 698c note that if classical boundary atom is 699c active the corresponding link atom is automatically 700c flagged as active in the first do loop because it 701c carries the same global index. The second do loop takes 702c care of the case when qm boundary atom is active using 703c link index constructed above 704c ------------------------------------------------ 705 nact1=0 706 do i=1,nfg 707 ii = int_mb(i_iqml+i-1) 708 do j=1,nact0 709 jj0 = int_mb(i_cons+j-1) 710 jj = int_mb(i_iqml+jj0-1) 711 aflag = .true. 712 if(aflag.and.(ii.eq.jj)) then 713 nact1 = nact1 + 1 714 int_mb(i_act+nact1-1) = i 715 goto 1 716 end if 717 end do 718 do j=1,nal 719 jj = int_mb(i_al+j-1) 720 if(ii.eq.jj) then 721 nact1 = nact1 + 1 722 int_mb(i_act+nact1-1) = i 723 goto 1 724 end if 725 end do 7261 continue 727 end do 728 729 if(nact1.ne.0) then 730 if (.not.rtdb_put(rtdb, 'geometry:actlist', 731 > mt_int,nact1,int_mb(i_act))) 732 > call errquit(pname,0, 733 > RTDB_ERR) 734 end if 735 736 if (.not.ma_free_heap(h_cons)) 737 $ call errquit(pname,h_cons, 738 & MA_ERR) 739 740 741 if (.not. ma_chop_stack(h_act) ) call errquit( 742 & pname//' ma_pop_stack ', 743 & 0, MA_ERR) 744 745 end 746 747 subroutine mm_restore_fixed_atoms(rtdb) 748 implicit none 749#include "rtdb.fh" 750#include "util.fh" 751#include "inp.fh" 752#include "mafdecls.fh" 753#include "errquit.fh" 754#include "nwc_const.fh" 755#include "geom.fh" 756 integer rtdb 757 758 character*30 pname 759 integer nact 760 integer i_cons,h_cons 761 integer ma_type 762 logical ignore 763 764 pname = "mm_restore_fixed_atoms" 765c write(*,*) pname 766 767 ignore = rtdb_delete(rtdb, 'geometry:actlist') 768 769 if (rtdb_ma_get(rtdb, 'qmmm:actlist', ma_type, 770 $ nact, h_cons)) then 771 if (.not. ma_get_index(h_cons, i_cons)) 772 $ call errquit(pname,h_cons, 773 & MA_ERR) 774 else 775 return 776 end if 777 778 if (.not.rtdb_delete(rtdb, 'qmmm:actlist')) 779 $ call errquit(pname,0, 780 & RTDB_ERR) 781 782 if (.not.rtdb_put(rtdb, 'geometry:actlist', 783 > mt_int,nact,int_mb(i_cons))) 784 $ call errquit(pname,0, 785 & RTDB_ERR) 786 787 if (.not.ma_free_heap(h_cons)) 788 $ call errquit(pname,h_cons, 789 & MA_ERR) 790 791 end 792