1 subroutine qmmm_bq_data_init(irtdb) 2 implicit none 3c 4#include "mafdecls.fh" 5#include "rtdb.fh" 6#include "bq.fh" 7#include "errquit.fh" 8#include "qmmm_bq_data.fh" 9 integer irtdb 10c 11 character*30 bq_expand 12 character*30 pname 13 14 pname = "qmmm_bq_data_init" 15 nbqs=0 16 nbqw=0 17 nbq=nbqs+nbqw 18 19 if (.not.rtdb_get(irtdb,'qmmm:bq_exclude',mt_int,1,bq_exclude)) 20 + call errquit(pname,bq_exclude, 21 & RTDB_ERR) 22 23 if (.not. rtdb_get(irtdb,'qmmm:bq_update',mt_int,1,bq_update)) 24 + bq_update = 0 25 26 readbq = .false. 27 if (rtdb_cget(irtdb,"qmmm:bq:load:file",1,bqfile_in)) then 28 if(bqfile_in.ne."none") then 29 readbq = .true. 30 if(bqfile_in.eq."default") then 31 call util_file_name("bq.ind",.false.,.false.,bqfile_in) 32 else 33 call util_file_name_resolve(bqfile_in,.false.) 34 end if 35 end if 36 end if 37 38 39 bqsave = .false. 40 if (rtdb_cget(irtdb,'qmmm:bq:save:file',1,bqfile_out)) then 41 if(bqfile_out.ne."none") then 42 bqsave = .true. 43 if(bqfile_out.eq."default") then 44 call util_file_name("bq.ind",.false.,.false.,bqfile_out) 45 else 46 call util_file_name_resolve(bqfile_out,.false.) 47 end if 48 end if 49 end if 50c 51 if (.not.rtdb_cget(irtdb,"qmmm:bq_expand",1,bq_expand)) 52 > bq_expand = "none" 53 if(bq_expand.eq."none") then 54 allbqs = .false. 55 allbqw = .false. 56 else if(bq_expand.eq."all") then 57 allbqs = .true. 58 allbqw = .true. 59 else if(bq_expand.eq."solute") then 60 allbqs = .true. 61 allbqw = .false. 62 else if(bq_expand.eq."solvent") then 63 allbqs = .false. 64 allbqw = .true. 65 else 66 call errquit(pname,"unknown option for charge expand",0,0) 67 end if 68c 69c if (.not.rtdb_get(irtdb,'qmmm:bq-solute',mt_log,1,allbqs)) 70c > allbqs = .false. 71cc 72c if (.not.rtdb_get(irtdb,'qmmm:bq-solvent',mt_log,1,allbqw)) 73c > allbqw = .false. 74 75 if(.not.bq_create("qmmm charges",bq_handle)) 76 + call errquit(pname//'Failed bq_create',0,CALC_ERR) 77 78 call qmmm_bq_data_load() 79c 80c force separate calculation of bq energy 81c --------------------------------------- 82 if (.not. rtdb_put(irtdb,'dft:bq_energy',mt_dbl,1,0.0d0)) 83 + call errquit(pname//'setting dft:bq_energy',0,CALC_ERR) 84c 85 86 end 87 88 subroutine qmmm_bq_data_load() 89 implicit none 90c 91#include "mafdecls.fh" 92#include "errquit.fh" 93#include "qmmm_bq_data.fh" 94#include "qmmm.fh" 95#include "qmmm_utils.fh" 96#include "qmmm_params.fh" 97#include "global.fh" 98#include "bq.fh" 99#include "util.fh" 100#include "rtdb.fh" 101 integer i,k 102 integer i_lb 103 integer i_mbq,h_mbq 104 integer nlink 105 integer il,j 106 double precision charge 107 integer mwa 108 character*32 pname 109 110 pname = "qmmm_bq_data_load" 111 112 if(qmmm_print_debug()) 113 > write(*,*) "in",pname 114 115 call qmmm_bq_data_dealloc() 116c 117c reset bq flags for solute and solvent 118c ------------------------------------- 119 if(allbqs) call mm_activate_bqszone() 120 if(allbqw) call mm_activate_bqwzone() 121c 122c exclude bq's as specified by qmmm input 123c --------------------------------------- 124 call mm_prune_bqzone() 125c 126c get numbers of solute and solvent bq's 127c -------------------------------------- 128 if(readbq) then 129 call qmmm_read_nbq() 130 else 131 call mm_get_tot_nbqs(nbqs) 132 call mm_get_tot_nbqw(nbqw) 133 end if 134 call mm_get_mwa(mwa) 135 nbq=nbqs+nbqw 136 mbqw = nbqw/mwa 137 138 if(nbq.eq.0) return 139 140c allocate bq structure 141c --------------------- 142 call qmmm_bq_data_alloc() 143 144c get global index for bq charges 145c ------------------------------- 146 if(readbq) then 147 call qmmm_read_ibq() 148 call mm_reset_solute_bqzone(nbqs,int_mb(i_ibq)) 149 call mm_reset_solvent_bqzone(nbqw,int_mb(i_ibq+nbqs)) 150 else 151 if(nbqs.ne.0) then 152 call mm_get_solute_bq_ind(nbqs,int_mb(i_ibq)) 153 end if 154 if(nbqw.ne.0) then 155 call mm_get_solvent_ind_bq(nbqw,int_mb(i_ibq+nbqs)) 156 end if 157 end if 158c 159c get solute bq coord and charges from MM 160c --------------------------------------- 161 if(nbqs.ne.0) then 162 call mm_get_solute_geom_bq(nbqs, 163 > int_mb(i_ibq), 164 > dbl_mb(i_cbq), 165 > dbl_mb(i_qbq)) 166 167 end if 168 169 if(qmmm_print_debug()) 170 > write(*,*) "finished setting up solute bq charges" 171c 172c append solvent bq coord and charges from MM 173c ------------------------------------------- 174 if(nbqw.ne.0) then 175 call mm_get_solvent_geom_bq(nbqw,int_mb(i_ibq+nbqs), 176 > dbl_mb(i_cbq+3*nbqs), 177 > dbl_mb(i_qbq+nbqs)) 178 end if 179 180 if(qmmm_print_debug()) 181 > write(*,*) "finished setting up solvent bq charges" 182 183c 184 if(qmmm_master()) then 185 write(*,*) "Total number of Bq charges",nbq 186 write(*,*) "Number of solute Bq charges",nbqs 187 write(*,*) "Number of solvent Bq charges",nbqw 188 end if 189c 190c save bq index into the file 191c --------------------------- 192 if(bqsave) 193 + call qmmm_save_bq_index() 194c 195c create and activate bq structure in bq module 196c --------------------------------------------- 197 if(.not.bq_pset(bq_handle,nbq,h_qbq,h_cbq)) 198 + call errquit(pname//'Failed bq_pset',0,CALC_ERR) 199 if(.not.bq_activate(bq_handle)) 200 + call errquit(pname//'Failed bq_activate',0,CALC_ERR) 201 charge = 0.0d0 202 do i=1,nbq 203 charge = charge + dbl_mb(i_qbq+i-1) 204 end do 205 if(ga_nodeid().eq.0) then 206 write(*,*) "Total Bq charge: ", charge 207 end if 208 209 if(util_print('bq_geom', print_high)) 210 + call bq_print_info(bq_handle) 211 212 if(ga_nodeid().eq.0) 213 > call qmmm_print_pdb_bq(nbq,"-bq.pdb", 214 > dbl_mb(i_cbq), 215 > dbl_mb(i_qbq)) 216 217 if(qmmm_print_debug()) then 218 call util_print_centered(6,"Bq information", 219 > 36, .true.) 220 write(*,'(/,A,T20,A,/)') "global index","coordinates(au)" 221 do i=1,nbqs+nbqw 222 write(*,'(I5,T20,3F12.6)') int_mb(i_ibq+i-1), 223 > (dbl_mb(i_cbq+3*(i-1)+k-1),k=1,3) 224 end do 225 end if 226 227 if(qmmm_print_debug()) 228 > write(*,*) "out",pname 229 230 end 231 232 subroutine qmmm_bq_data_reload() 233 implicit none 234c 235#include "mafdecls.fh" 236#include "errquit.fh" 237#include "qmmm_bq_data.fh" 238#include "qmmm.fh" 239 integer counter 240 character*32 pname 241 save counter 242 DATA counter /0/ 243 244 pname = "qmmm_bq_data_reload" 245 246 if(qmmm_print_debug()) 247 > write(*,*) "in",pname 248 249 counter = counter + 1 250 if(counter.eq.bq_update) then 251 call qmmm_bq_data_load() 252 counter = 0 253 end if 254 255 if(qmmm_print_debug()) 256 > write(*,*) "out",pname 257 258 return 259 end 260 261 subroutine qmmm_bq_data_update_active() 262 implicit none 263c 264#include "mafdecls.fh" 265#include "errquit.fh" 266#include "qmmm_bq_data.fh" 267#include "qmmm.fh" 268#include "qmmm_utils.fh" 269#include "qmmm_params.fh" 270#include "global.fh" 271#include "bq.fh" 272#include "util.fh" 273 integer i,k 274 character*32 pname 275 276 pname = "qmmm_bq_data_create_active" 277 278 if(qmmm_print_debug()) 279 > write(*,*) "in",pname 280c 281 do i=1,nbq 282 log_mb(i_abq+i-1)= .false. 283 end do 284c 285c get active list of bqs 286c ------------------------ 287 if(nbqs.ne.0) 288 > call qmmm_cons_get_actmaps(nbqs,int_mb(i_ibq), 289 > log_mb(i_abq)) 290 291 if(nbqw.ne.0) 292 > call qmmm_cons_get_actmapw(nbqw,int_mb(i_ibq+nbqs), 293 > log_mb(i_abq+nbqs)) 294 295 nbqa = 0 296 do i=1,nbq 297 if(log_mb(i_abq+i-1)) then 298 nbqa = nbqa + 1 299 end if 300 end do 301 302 if(ga_nodeid().eq.0) then 303 write(*,*) "Total number of active Bq charges ",nbqa 304 end if 305 306 if(qmmm_print_debug()) 307 > write(*,*) "out",pname 308 309 end 310 311 subroutine qmmm_bq_data_reset_active() 312 implicit none 313c 314#include "mafdecls.fh" 315#include "errquit.fh" 316#include "qmmm_bq_data.fh" 317#include "qmmm.fh" 318#include "qmmm_utils.fh" 319#include "qmmm_params.fh" 320#include "global.fh" 321#include "bq.fh" 322#include "util.fh" 323 integer i,k 324 character*32 pname 325 326 pname = "qmmm_bq_data_create_active" 327 328 if(qmmm_print_debug()) 329 > write(*,*) "in",pname 330c 331 do i=1,nbq 332 log_mb(i_abq+i-1)= .true. 333 end do 334 nbqa = nbq 335 do i=1,nbq 336 if(log_mb(i_abq+i-1)) then 337 nbqa = nbqa + 1 338 end if 339 end do 340 341 if(ga_nodeid().eq.0) then 342 write(*,*) "number of active bqs ",nbqa 343 end if 344 345 if(qmmm_print_debug()) 346 > write(*,*) "out",pname 347 348 end 349 350 subroutine qmmm_bq_print_info() 351 implicit none 352c 353#include "mafdecls.fh" 354#include "errquit.fh" 355#include "qmmm_bq_data.fh" 356#include "qmmm.fh" 357#include "qmmm_utils.fh" 358#include "qmmm_params.fh" 359#include "global.fh" 360#include "bq.fh" 361#include "util.fh" 362 integer i,k 363 character*32 pname 364 365 pname = "qmmm_bq_print_info" 366 367 if(qmmm_print_debug()) 368 > write(*,*) "in",pname 369 370 call util_print_centered(6,"Bq information", 371 > 36, .true.) 372 write(*,'(/,A,T20,A,/)') "global index","coordinates(au)" 373 do i=1,nbqs+nbqw 374 write(*,'(I5,T20,3F12.6)') int_mb(i_ibq+i-1), 375 > (dbl_mb(i_cbq+3*(i-1)+k-1),k=1,3) 376 end do 377 378 if(qmmm_print_debug()) 379 > write(*,*) "out",pname 380 381 end 382 383 subroutine qmmm_bq_coord_update() 384 implicit none 385c 386#include "mafdecls.fh" 387#include "errquit.fh" 388#include "qmmm_bq_data.fh" 389#include "qmmm.fh" 390#include "global.fh" 391c 392 integer i 393 if(nbqs.ne.0) then 394 call mm_get_solute_coord(nbqs, 395 > int_mb(i_ibq), 396 > dbl_mb(i_cbq)) 397 398 end if 399 400c 401c append solvent bq coord and charges from MM 402c ------------------------------------------- 403 if(nbqw.ne.0) then 404 call mm_get_solvent_coord(nbqw,int_mb(i_ibq+nbqs), 405 > dbl_mb(i_cbq+3*nbqs)) 406 end if 407 408 if(qmmm_print_debug()) 409 > write(*,*) "finshed updating up bq coords" 410 411 end 412 413 subroutine qmmm_bq_data_alloc() 414 implicit none 415c 416#include "mafdecls.fh" 417#include "errquit.fh" 418#include "qmmm_bq_data.fh" 419 integer i 420 if(nbq.eq.0) return 421 422 if(.not.ma_alloc_get(mt_log,nbq,'abq',h_abq,i_abq)) 423 + call errquit('qmmm: Failed to allocate memory for abq',nbq, 424 + MA_ERR) 425 do i=1,nbq 426 log_mb(i_abq+i-1)=.true. 427 end do 428 429 if(.not.ma_alloc_get(mt_dbl,3*nbq,'cbq',h_cbq,i_cbq)) 430 + call errquit('qmmm: Failed to allocate memory for cbq',3*nbq, 431 + MA_ERR) 432 call dfill(3*nbq,0.d0,dbl_mb(i_cbq),1) 433 if(.not.ma_alloc_get(mt_dbl,nbq,'qbq',h_qbq,i_qbq)) 434 + call errquit('qmmm: Failed to allocate memory for qbq',nbq, 435 & MA_ERR) 436 call dfill(nbq,0.d0,dbl_mb(i_qbq),1) 437 if(.not.ma_alloc_get(mt_int,nbq,'ibq',h_ibq,i_ibq)) 438 + call errquit('qmmm: Failed to allocate memory for ibq',nbq, 439 & MA_ERR) 440 call ifill(nbq,0,int_mb(i_ibq),1) 441 442 if(.not.ma_alloc_get(mt_int,nbq,'i_ibqa',h_ibqa,i_ibqa)) 443 + call errquit('qmmm: Failed to allocate memory for i_ibqa',nbq, 444 & MA_ERR) 445 call ifill(nbq,0,int_mb(i_ibqa),1) 446 447 end 448 449 subroutine qmmm_bq_data_dealloc() 450 implicit none 451c 452#include "mafdecls.fh" 453#include "errquit.fh" 454#include "qmmm_bq_data.fh" 455#include "bq.fh" 456 logical ignore 457 458 if(nbq.eq.0) return 459 460 if(.not.ma_free_heap(h_ibqa)) 461 + call errquit('qmmm:memory deallocation for ibqa',nbqa, 462 & MA_ERR) 463 464 if(.not.ma_free_heap(h_ibq)) 465 + call errquit('qmmm:memory deallocation for ibq',nbq, 466 & MA_ERR) 467 468 if(.not.ma_free_heap(h_qbq)) 469 + call errquit('qmmm:memory deallocation for qbq',nbq, 470 & MA_ERR) 471 472 if(.not.ma_free_heap(h_cbq)) 473 + call errquit('qmmm: memory deallocation for cbq',3*nbq, 474 + MA_ERR) 475 476 if(.not.ma_free_heap(h_abq)) 477 + call errquit('qmmm: memory deallocation for cbq',3*nbq, 478 + MA_ERR) 479 480 nbqs=0 481 nbqw=0 482 nbq=0 483 484 ignore = bq_destroy(bq_handle) 485 486 487 end 488 489 function qmmm_get_nbq() 490 implicit none 491#include "qmmm_bq_data.fh" 492 493 integer qmmm_get_nbq 494 495 qmmm_get_nbq = nbqs+nbqw 496 497 end 498 499 function qmmm_get_nbqa() 500 implicit none 501#include "qmmm_bq_data.fh" 502 503 integer qmmm_get_nbqa 504 505 qmmm_get_nbqa = nbqa 506 507 end 508 509 function qmmm_get_nbqs() 510 implicit none 511#include "qmmm_bq_data.fh" 512 513 integer qmmm_get_nbqs 514 515 qmmm_get_nbqs = nbqs 516 517 end 518 519 function qmmm_get_nbqw() 520 implicit none 521#include "qmmm_bq_data.fh" 522 523 integer qmmm_get_nbqw 524 525 qmmm_get_nbqw = nbqw 526 527 end 528 529 function qmmm_get_h_qbq() 530 implicit none 531#include "qmmm_bq_data.fh" 532 integer qmmm_get_h_qbq 533 534 qmmm_get_h_qbq = h_qbq 535 536 end 537 538 function qmmm_get_i_qbq() 539 implicit none 540#include "qmmm_bq_data.fh" 541 integer qmmm_get_i_qbq 542 543 qmmm_get_i_qbq = i_qbq 544 545 end 546 547 function qmmm_get_h_cbq() 548 implicit none 549#include "qmmm_bq_data.fh" 550 integer qmmm_get_h_cbq 551 552 qmmm_get_h_cbq = h_cbq 553 554 end 555 556 function qmmm_get_i_cbq() 557 implicit none 558#include "qmmm_bq_data.fh" 559 integer qmmm_get_i_cbq 560 561 qmmm_get_i_cbq = i_cbq 562 563 end 564 565 function qmmm_get_i_ibq() 566 implicit none 567#include "qmmm_bq_data.fh" 568 integer qmmm_get_i_ibq 569 570 qmmm_get_i_ibq = i_ibq 571 572 end 573 574 function qmmm_get_i_ibqa() 575 implicit none 576#include "qmmm_bq_data.fh" 577 integer qmmm_get_i_ibqa 578 579 qmmm_get_i_ibqa = i_ibqa 580 581 end 582 583 function qmmm_get_h_cbqs() 584 implicit none 585#include "qmmm_bq_data.fh" 586 integer qmmm_get_h_cbqs 587 588 qmmm_get_h_cbqs = h_cbqs 589 590 end 591 592 function qmmm_get_i_cbqs() 593 implicit none 594#include "qmmm_bq_data.fh" 595 integer qmmm_get_i_cbqs 596 597 qmmm_get_i_cbqs = i_cbqs 598 599 end 600 601 function qmmm_get_h_cbqw() 602 implicit none 603#include "qmmm_bq_data.fh" 604 integer qmmm_get_h_cbqw 605 606 qmmm_get_h_cbqw = h_cbqw 607 608 end 609 610 function qmmm_get_i_cbqw() 611 implicit none 612#include "qmmm_bq_data.fh" 613 integer qmmm_get_i_cbqw 614 615 qmmm_get_i_cbqw = i_cbqw 616 617 end 618 619 function qmmm_get_i_qbqw() 620 implicit none 621#include "qmmm_bq_data.fh" 622 integer qmmm_get_i_qbqw 623 624 qmmm_get_i_qbqw = i_qbqw 625 626 end 627 628 function qmmm_get_bq_exclude() 629 implicit none 630#include "qmmm_bq_data.fh" 631 632 integer qmmm_get_bq_exclude 633 634 qmmm_get_bq_exclude = bq_exclude 635 636 end 637 638 subroutine qmmm_read_ibq() 639 implicit none 640#include "stdio.fh" 641#include "errquit.fh" 642#include "qmmm.fh" 643#include "util.fh" 644#include "inp.fh" 645#include "mafdecls.fh" 646#include "global.fh" 647#include "msgids.fh" 648#include "qmmm_bq_data.fh" 649 650c local variables 651 integer nf 652 character*30 pname 653 654 integer i 655 integer fn 656 integer anbq,anbqs,anbqw 657 658 pname = "qmmm_read_ibq" 659 660 if(qmmm_master()) then 661 if(.not.qmmm_get_io_unit(fn)) 662 > call errquit("cannot get file number",0,0) 663 open(unit=fn,status="old",form="formatted",file=bqfile_in) 664 read(fn,*) anbqs,anbqw,anbq 665 do i=1,anbq 666 read(fn,*) int_mb(i_ibq+i-1) 667 end do 668 close(fn) 669 call util_print_centered(luout, 670 > "Bq index has been read from external file "// 671 > bqfile_in, 672 > 36, .true.) 673 end if 674 call ga_brdcst(msg_qmmm_ind,int_mb(i_ibq), 675 > nbq*ma_sizeof(mt_int,1,mt_byte),0) 676 677 call ga_sync() 678 679 end 680 681 subroutine qmmm_read_nbq() 682 implicit none 683#include "errquit.fh" 684#include "qmmm.fh" 685#include "util.fh" 686#include "inp.fh" 687#include "mafdecls.fh" 688#include "global.fh" 689#include "msgids.fh" 690#include "qmmm_bq_data.fh" 691 692c local variables 693 integer nf 694 695 integer i 696 integer fn 697 698 if(qmmm_master()) then 699 if(.not.qmmm_get_io_unit(fn)) 700 > call errquit("cannot get file number",0,0) 701 open(unit=fn,status="old",form="formatted",file=bqfile_in) 702 read(fn,*) nbqs,nbqw 703 close(fn) 704 end if 705 706 call ga_brdcst(msg_qmmm_nbqs,nbqs,ma_sizeof(mt_int,1,mt_byte),0) 707 call ga_brdcst(msg_qmmm_nbqw,nbqw,ma_sizeof(mt_int,1,mt_byte),0) 708 709 710 end 711 712 subroutine qmmm_save_bq_index() 713 implicit none 714#include "stdio.fh" 715#include "errquit.fh" 716#include "qmmm.fh" 717#include "util.fh" 718#include "global.fh" 719#include "inp.fh" 720#include "qmmm_bq_data.fh" 721#include "mafdecls.fh" 722 723c local variables 724 integer nf 725 character*(nw_max_path_len) filename 726 727 integer i 728 integer fn 729 730 if(qmmm_master()) then 731 732 if(.not.qmmm_get_io_unit(fn)) 733 > call errquit("cannot get file number",0,0) 734 735 open(unit=fn,status="unknown",form="formatted",file=bqfile_out) 736 737 write(fn,*) nbqs,nbqw,nbq 738 do i=1,nbq 739 write(fn,*) int_mb(i_ibq+i-1) 740 end do 741 742 call util_flush(fn) 743 close(fn) 744 745 call util_print_centered(luout, 746 > "Bq index has been saved to external file "// 747 > bqfile_out, 748 > 36, .true.) 749 750 end if 751 752 call ga_sync() 753 754 end 755 756 subroutine qmmm_bq_push_grad(nt, 757 > g) 758 implicit none 759 760#include "qmmm_bq_data.fh" 761#include "qmmm_params.fh" 762#include "global.fh" 763#include "mafdecls.fh" 764#include "errquit.fh" 765 766 integer nt 767 double precision g(3*nt) 768c 769 integer i_ibqg,h_ibqg 770 integer i,j,k,ia 771 integer ii,jj 772 integer nbqsa,nbqwa,mwa 773 character*32 pname 774 775 pname = "qmmm_bq_push_grad" 776 if(.not.ma_push_get(mt_int,nbqa,'ibqg',h_ibqg,i_ibqg)) 777 + call errquit( pname//'Failed to allocate memory for ibqs', 778 + nbqa, MA_ERR) 779 780 781 nbqsa = 0 782 do i=1,nbqs 783 if(log_mb(i_abq+i-1)) then 784 nbqsa = nbqsa + 1 785 int_mb(i_ibqg+nbqsa-1) = int_mb(i_ibq+i-1) 786 end if 787 end do 788 789 call mm_get_mwa(mwa) 790 nbqwa = 0 791 do i=1,nbq-nbqs,mwa 792 if(log_mb(i_abq+nbqs+i-1)) then 793 nbqwa = nbqwa + 1 794 int_mb(i_ibqg+nbqsa+nbqwa-1) = int_mb(i_ibq+nbqs+(i-1)/mwa) 795 end if 796 end do 797 798 799 if(nbqsa.ne.0) 800 > call mm_add_solute_force(nbqsa, 801 > int_mb(i_ibqg), 802 > g) 803 804 805 if(nbqwa.ne.0) 806 > call mm_add_solvent_force(mwa*nbqwa, 807 > int_mb(i_ibqg+nbqsa), 808 > g(3*nbqsa+1)) 809 810 811 if(.not.ma_pop_stack(h_ibqg)) 812 + call errquit( pname//'Failed to deallocate memory for ibqg', 813 + nbqa, MA_ERR) 814 815 816 end 817 818 subroutine qmmm_get_cqbqa(nt, 819 > c, 820 > q) 821 implicit none 822 823#include "qmmm_bq_data.fh" 824#include "global.fh" 825#include "mafdecls.fh" 826#include "errquit.fh" 827#include "qmmm_params.fh" 828 829 integer nt 830 double precision c(3,nt) 831 double precision q(nt) 832c 833 integer i,j,k,ia 834 integer ii,jj 835 836 837 ia = 0 838 do i=1,nbqs 839 if(log_mb(i_abq+i-1)) then 840 ia = ia + 1 841 q(ia) = dbl_mb(i_qbq+i-1) 842 do k=1,3 843 c(k,ia) = dbl_mb(i_cbq+(i-1)*3+k-1) 844 end do 845 end if 846 end do 847 848 do i=nbqs+1,nbq 849 if(log_mb(i_abq+i-1)) then 850 ia = ia + 1 851 q(ia) = dbl_mb(i_qbq+i-1) 852 do k=1,3 853 c(k,ia) = dbl_mb(i_cbq+(i-1)*3+k-1) 854 end do 855 end if 856 end do 857 858 859 end 860 861c $Id$ 862