1 subroutine bsse_input(rtdb) 2 implicit none 3c 4c This subroutine read the input lines in bsse 5c to get info about monomers that build the 6c super molecule. nmon identify monomers, mon(nmon), 7c each monomer have mon_atm(nmon). Name of monomers is 8c allocated in mon_name(nmon) 9c 10#include "stdio.fh" 11#include "errquit.fh" 12#include "mafdecls.fh" 13#include "inp.fh" 14#include "rtdb.fh" 15#include "geom.fh" 16c here are defined common variables 17#include "bsse_common.fh" 18c 19 integer rtdb ! [input] 20c 21 integer geom 22 integer atm_tot ! to check total atoms 23 integer qtot ! to print total charge 24 integer nopen ! to print spin multiplicity 25 integer i, j,l,k 26 integer nfield 27 integer nopt 28 parameter (nopt = 9) 29 integer ind 30c 31 character*80 buf 32 character*255 field 33 character*18 opt (nopt) 34c 35 logical status 36 logical do_bsse 37c 38c: 39 logical bsse_rtdb_store 40 external bsse_rtdb_store 41 logical bsse_rtdb_load 42 external bsse_rtdb_load 43c 44 data opt /'end', 'on', 'off','tidy','charge', 'input', 45 & 'input_wghost', 'mon', 'mult'/ 46 47c 48c ------------------welcome------------------------ 49c 50 buf = ' ' 51 write(buf,*) ' Input BSSE Module - Counter Poise Approach' 52 write(LuOut,*) 53 write(LuOut,*) 54 call util_print_centered(LuOut,buf,40,.true.) 55 write(LuOut,*) 56c 57c ------------------------------------------------- 58c ----- get info supermolecule geometry ------- 59c ------------------------------------------------- 60c 61 if (.not. rtdb_cget(rtdb,'geometry', 1, spr_name)) 62 $ spr_name = 'geometry' 63c 64 65 if (.not. geom_create(geom, spr_name)) 66 $ call errquit('bsse_input: geom_create failed !', 0,GEOM_ERR) 67c 68 69 if (.not. geom_rtdb_load(rtdb, geom, spr_name)) 70 $ call errquit('bsse_input: no geometry load form rtdb', 0, 71 $ GEOM_ERR) 72 73 if (.not. geom_ncent(geom, natoms)) call errquit 74 $ ('bsse_input: geom_ncent ?', 0, GEOM_ERR) 75c ------------------------------------------------- 76c ----- reading input fields -------- 77c ------------------------------------------------- 78c 79 80 call inp_set_field(0) ! goto the begin of line 81 82c 83c: preliminaries 84c 85 qtot = 0 86 atm_tot = 0 87 nmon = 0 88 do_bsse = .true. 89 90c 91 92 call ifill(mx_atm,1,mmon,1) ! multiplicity default 93 call dfill(mx_atm,0.0d0,qmon,1) ! charge default 94 call cfill(mx_atm,' ',input,1) ! charge default 95c 96 100 if(.not.inp_read()) 97 $ call errquit('bsse_input: unexpected eof ',911,INPUT_ERR ) 98c 99 nfield = inp_n_field() 100c 101 150 if (.not.inp_a(field)) 102 $ call errquit('bsse_input: failed to read field',911,INPUT_ERR ) 103c 104 if (.not. inp_match(nopt, .false., field, opt, ind)) 105 $ goto 10 106 107 goto (900, 850, 800, 700, 600, 500, 400, 300, 200) ind 108 109c 110c: none 111c 112 113 10 write(LuOut,20) 114 115 20 format 116 117 $(/' valid bsse structure input : '/ 118 $ ' bsse '/ 119 $ ' mon <character name fragment> <integer list atoms>'/ 120 $ ' input <input line> '/ 121 $ ' input_wghost <input line> '/ 122 $ ' charge <double charge> '/ 123 $ ' mult <integer multiplicity> '/ 124 $ ' off '/ 125 $ ' on '/ 126 $ ' end '/ 127 $/) 128c 129 call errquit('bsse_input: unknown directive', 911,INPUT_ERR) 130c 131c 132c: mon 133c 134 135 300 continue 136c 137 nmon = nmon + 1 138c 139 mon_name(nmon) = field ! took the first field however use to name 140c 141 if (.not.inp_a(mon_name(nmon))) 142 $ call errquit 143 $ ('bsse_input: failed to read name field',911,INPUT_ERR) 144c 145c Read the atom numbers and count the number of atoms as we go along. 146c If we read something else than an integer it might be the next 147c keyword on the line. So, leave the loop, and check if we have reached 148c the end of the line. If we are not at the end of the line goto 150 149c to read the next keyword, otherwise goto 100 to read the next line. 150c 151 i = 0 152 do while (inp_i(mon(nmon,i+1))) 153 i = i + 1 154 mon_atm(nmon) = i 155 atm_tot = atm_tot + 1 156 enddo 157 if (inp_cur_field().lt.nfield) goto 150 158c 159 go to 100 160c 161c: mult 162c 163 200 continue 164 if(nmon.eq.0) goto 10 165 if (.not. inp_i(mmon(nmon) )) call errquit 166 $ ('bsse_input: failed reading monomer multiplicity', 167 $ nmon,INPUT_ERR) 168 if (mmon(nmon).eq.0) call errquit 169 $ ('bsse_input: invalid multiplicity ',mmon(nmon), 170 $ INPUT_ERR) 171 if (inp_cur_field().lt.nfield) goto 150 172 goto 100 173c 174c: input_wghost 175c 176 400 continue 177c 178 if(nmon.eq.0) goto 10 179 180 i=(nmon)*2 181 182 if (.not. inp_a(input(i) )) call errquit 183 $ ('bsse_input: failed reading input [input]',911,INPUT_ERR) 184 if (inp_cur_field().lt.nfield) goto 150 185 186 go to 100 187c 188c: input 189c 190 500 continue 191 if(nmon.eq.0) goto 10 192 193 i=(nmon-1)*2+1 194 195 if (.not. inp_a(input(i))) 196 $ call errquit 197 $ ('bsse_input: failed reading input [input]',911,INPUT_ERR) 198c 199 if (inp_cur_field().lt.nfield) goto 150 200c 201 go to 100 202c 203 204 600 continue 205c 206c: charge 207c 208 if(nmon.eq.0) goto 10 209c 210 if(.not. inp_f( qmon(nmon))) 211 $ call errquit('bsse_input: reading monomer charge',911,INPUT_ERR) 212c 213 if (inp_cur_field().lt.nfield) goto 150 214 215 go to 100 216c 217c: tidy 218c clean database of any bsse info 219 700 continue 220c 221c: off 222c 223 if(.not. rtdb_delete(rtdb, 'bsse')) 224 $ call errquit('bsse_input: cannot clean database',911,RTDB_ERR) 225 goto 100 226c 227 800 continue 228 229 do_bsse=.false. 230 buf = ' ' 231 write(buf,*) ' Any BSSE operations are off' 232 write(LuOut,*) 233 write(LuOut,*) 234 call util_print_centered(LuOut,buf,40,.true.) 235 write(LuOut,*) 236 237 goto 100 238c 239c 240 850 continue 241 242 do_bsse=.true. 243 if(.not.bsse_rtdb_load(rtdb)) 244 $ call errquit('bsse_input: load data input in db',911,RTDB_ERR) 245c 246 atm_tot= natoms 247c 248 goto 100 249c 250c: end 251c 252 900 continue 253 254 if(do_bsse) then 255 256 goto 1000 257 258 else 259 260 goto 1100 261 262 endif 263c 264c: done 265c 266 1000 continue 267 268c ------------------------------------------------------ 269c check : total atoms 270c ------------------------------------------------------ 271 272 if (atm_tot.ne.natoms) 273c $ goto 10 274 $ call errquit 275 $ ('bsse_input: number of atoms is wrong',911,INPUT_ERR) 276 277c ------------------------------------------------------ 278c check : dont repeat atoms 279c ------------------------------------------------------ 280 281 do j = 1, nmon 282 283 do i = 1, mon_atm(j) 284 285 do l = j, nmon 286 287 do k = 1, mon_atm(l) 288 289 if((i.eq.k).and.(j.eq.l)) then 290 291 status=.true. 292 293 elseif (mon(j,i).eq.mon(l,k)) then 294 295 call errquit 296 $ ('bsse_input: there are some atoms repeated',911,INPUT_ERR) 297 endif 298 enddo 299 enddo 300 enddo 301 enddo 302 303c---------------------------------------------------- 304c check the atoms are correct 305c---------------------------------------------------- 306 do j = 1, nmon 307 do i = 1, mon_atm(j) 308 309 if( mon(j,i).gt.natoms) then 310 call errquit 311 $ ('bsse_input: incorrect atom number',mon(j,i),INPUT_ERR) 312 endif 313 enddo 314 enddo 315 316c ------------------------------------------------------ 317c check : total charge 318c ------------------------------------------------------ 319 320 do j = 1, nmon 321 322 qtot = qtot + qmon(j) 323 324 enddo 325c 326c ------------------------------------------------------ 327c check : unpaired electrons 328c ------------------------------------------------------ 329c 330 nopen = 0 331 do j = 1, nmon 332 nopen = nopen + abs(mmon(j)) - 1 333 enddo 334c 335 write(LuOut, 60) spr_name, nmon, 336 $ mod(nopen,2)+1, nopen+1, qtot 337 60 format(/ 338 $ ' supermolecule geometry name = ', a50/ 339 $ ' number of monomers = ', i4/ 340 $ ' total multiplicity = ', i4,' to ',i4/ 341 $ ' total charge = ', i4/) 342c 343 344 write(LuOut, *) 345 $ ' atoms for each monomer ' 346 347 do j=1, nmon 348 349 write(LuOut, 70) mon_name(j) 350 70 format(/ 351 $ ' monomer ', a10,' : ',$) 352c 353 354 do i=1, mon_atm(j) 355 356 write(LuOut, 90) mon(j,i) 357 358 90 format (i3,$) 359 360 enddo 361 enddo 362 363 write(LuOut,*) 364 call util_flush(luout) 365 366c 367 go to 1100 368c 369c: db is done and geometry is destroyed 370c 371c 1100 continue 372c 373 374 1100 continue 375c 376 if (do_bsse) then 377c 378 if (.not. rtdb_put(rtdb,'bsse',mt_log,1,.true.)) 379 $ call errquit('bsse_input: rtdb_put failed',911,RTDB_ERR) 380c 381 else 382c 383 if (.not. rtdb_put(rtdb,'bsse',mt_log,1,.false.)) 384 $ call errquit('bsse_input: rtdb_put failed',911,RTDB_ERR) 385c 386 endif 387c 388 if(.not. geom_destroy(geom)) 389 $ call errquit('bsse_input: geom_destroy failed', 911,RTDB_ERR) 390c 391 if(.not.bsse_rtdb_store(rtdb)) 392 $ call errquit('bsse_input: store data input in db',911,RTDB_ERR) 393c 394 return 395c 396 end 397C> 398C> \brief Initialize monomer calculation 399C> 400C> Initialize the monomer calculation by modifying the contents of 401C> the RunTime Data Base for the current calculation. 402C> 403 subroutine bsse_param(rtdb, mult, charge, j_mon_name, 404 & i_input,theory) 405 implicit none 406 integer rtdb !< [Input] The RTDB handle 407 character*(*) j_mon_name !< [Input] Monomer name 408 character*(*) i_input !< [Input] Line of input for monomer calculation 409 character*(*) theory !< [Input] The theory to apply 410 integer mult !< [Input] Monomer spin multiplicity 411 double precision charge !< [Input] Monomer charge 412 logical first_j 413 character*255 vec_dbi, vec_dbo,tmp 414 integer lentheo, lenname 415#include "rtdb.fh" 416#include "mafdecls.fh" 417#include "errquit.fh" 418#include "inp.fh" 419#include "stdio.fh" 420c 421 logical nw_inp_from_character 422 external nw_inp_from_character 423c 424 if (.not. rtdb_get(rtdb,'bsse:first_j',mt_log,1,first_j)) 425 $ first_j=.true. 426c 427c if (.not. rtdb_cget(rtdb,'theory', 1, theory)) 428c & call errquit('bsse_param: get theory',0) 429c 430 lenname = inp_strlen(j_mon_name) 431c 432 lentheo = inp_strlen(theory) 433c 434c: multiplicity 435c: density methods 436 if ( theory(1:lentheo).eq.'dft' .or. 437 $ theory(1:lentheo).eq.'tddft') then 438 if (.not. rtdb_put(rtdb, 'dft:mult', mt_int, 1, mult)) 439 $ call errquit('bsse_param: rtdb_put of mult failed', 440 $ 0,RTDB_ERR ) 441c: wavefuntion methods 442 443 elseif( theory(1:lentheo).ne.'dft' .and. 444 $ theory(1:lentheo).ne.'tddft') then 445 if (.not. rtdb_put(rtdb,'scf:nopen', MT_INT, 1, mult-1)) 446 $ call errquit('bsse_param: rtdb_put of nopen failed', 447 $ 0,RTDB_ERR) 448 endif 449c 450 if (charge .ne. -999.0d0) then 451c 452 if (.not. rtdb_put(rtdb,'charge',mt_dbl,1,charge)) 453 $ call errquit('bsse_param:setting charge?',911, RTDB_ERR) 454c 455 end if 456c 457 if (j_mon_name .ne. ' ') then 458c 459 if (.not. rtdb_cput(rtdb,'geometry',1,j_mon_name)) 460 $ call errquit('bsse_param: setting geometry?',911,RTDB_ERR) 461c 462 end if 463c 464 if (first_j) then 465c 466 tmp = 'atomic' 467c 468 else 469c 470 tmp = ' ' 471 tmp = j_mon_name(1:lenname)//'.bsse.movecs' 472c 473 endif 474c MP2 SCF DFT vectors 475c 476 vec_dbi = ' ' 477 vec_dbo = ' ' 478c 479c if(theory(1:lentheo).ne.'scf' .or. theory(1:lentheo).ne.'dft' 480c & .or. theory(1:lentheo).ne.'mcscf') then 481 if(theory(1:lentheo).eq.'dft') then 482 vec_dbi= theory(1:lentheo)//':input vectors' 483 vec_dbo= theory(1:lentheo)//':output vectors' 484c 485 elseif (theory(1:lentheo).eq.'mcscf') then 486 vec_dbi= theory(1:lentheo)//':input vectors' 487 vec_dbo= theory(1:lentheo)//':output vectors' 488c 489 else 490 vec_dbi = 'scf:input vectors' 491 vec_dbo = 'scf:output vectors' 492 493 endif 494c 495c 496 if (.not. rtdb_cput(rtdb, vec_dbi, 1, tmp)) 497 & call errquit('bsse_param: input_vectors',0,RTDB_ERR) 498c 499 tmp = ' ' 500 tmp = j_mon_name(1:lenname)//'.bsse.movecs' 501c 502 if (.not. rtdb_cput(rtdb, vec_dbo, 1,tmp)) 503 & call errquit('bsse_param: output_vectors',0,RTDB_ERR) 504c 505c 506 if (i_input .ne. ' ') then 507 508 if (.not. nw_inp_from_character(rtdb,i_input)) 509 $ call errquit('bsse_param: error processing input string', 510 & 060,INPUT_ERR) 511 512 endif 513c 514 end 515c 516c 517 logical function bsse_rtdb_store(rtdb) 518 implicit none 519#include "rtdb.fh" 520#include "errquit.fh" 521c#include "tcgmsg.fh" 522#include "bsse_common.fh" 523#include "mafdecls.fh" 524#include "stdio.fh" 525 526c the propose is store bsse's sets into db 527 528 integer rtdb ![input] 529c integer itmp ,k 530 character*255 ctmp 531c 532 ctmp = 'bsse:natoms' 533 if (.not. rtdb_put( rtdb, ctmp, mt_int, 1,natoms )) 534 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 535c 536c 537 ctmp = 'bsse:nmon' 538 if (.not. rtdb_put( rtdb, ctmp, mt_int, 1,nmon )) 539 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 540c 541 ctmp = 'bsse:mon_name' 542 if(.not.rtdb_cput(rtdb, ctmp, nmon, mon_name )) 543 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 544c 545 ctmp = 'bsse:spr_name' 546 if(.not.rtdb_cput(rtdb, ctmp, 1, spr_name )) 547 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 548c 549 ctmp = 'bsse:mon_atm' 550 if (.not. rtdb_put( rtdb, ctmp, mt_int, nmon, mon_atm )) 551 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 552c 553 ctmp = 'bsse:mon' 554 if (.not. rtdb_put( rtdb, ctmp, mt_int,mx_atm*mx_atm, mon)) 555 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 556c 557 ctmp = 'bsse:qmon' 558 if (.not. rtdb_put( rtdb, ctmp, mt_dbl ,nmon, qmon)) 559 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 560c 561 ctmp = 'bsse:mmon' 562 if (.not. rtdb_put( rtdb, ctmp, mt_int ,nmon, mmon)) 563 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 564c 565 ctmp= 'bsse:input' 566 if(.not.rtdb_cput(rtdb, ctmp, nmon*2, input)) 567 $ call errquit('bsse_rtdb_store: rtdb_put failed',0,RTDB_ERR) 568c 569 bsse_rtdb_store = .true. 570 return 571 end 572c 573 logical function bsse_rtdb_load(rtdb) 574 implicit none 575#include "rtdb.fh" 576#include "errquit.fh" 577#include "stdio.fh" 578c#include "tcgmsg.fh" 579#include "bsse_common.fh" 580#include "mafdecls.fh" 581 integer rtdb ![input] 582 character*255 ctmp 583c 584c the propose is load bsse's sets from db 585c 586 ctmp = 'bsse:natoms' 587 if (.not.rtdb_get( rtdb, ctmp, mt_int, 1, natoms)) 588 $ call errquit('bsse_rtdb_load: rtdb_put failed',0,RTDB_ERR) 589c 590 ctmp = 'bsse:nmon' 591 if(.not.rtdb_get( rtdb, ctmp, mt_int, 1 , nmon)) 592 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 593c 594 ctmp = 'bsse:mon_name' 595 if(.not.rtdb_cget( rtdb, ctmp, nmon, mon_name)) 596 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 597c 598 ctmp = 'bsse:spr_name' 599 if(.not.rtdb_cget( rtdb, ctmp, 1, spr_name)) 600 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 601c 602 ctmp = 'bsse:mon_atm' 603 if(.not.rtdb_get( rtdb, ctmp, mt_int, nmon,mon_atm)) 604 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 605c 606 ctmp = 'bsse:mon' 607 if (.not.rtdb_get( rtdb, ctmp, mt_int, mx_atm*mx_atm,mon)) 608 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 609c 610 ctmp = 'bsse:qmon' 611 if (.not.rtdb_get( rtdb, ctmp, mt_dbl, nmon, qmon )) 612 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 613c 614 ctmp = 'bsse:mmon' 615 if (.not. rtdb_get( rtdb, ctmp, mt_int ,nmon, mmon)) 616 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 617c 618 ctmp = 'bsse:input' 619 if(.not.rtdb_cget( rtdb, ctmp, nmon*2, input)) 620 $ call errquit('bsse_rtdb_load: rtdb_get failed',0,RTDB_ERR) 621c 622 bsse_rtdb_load = .true. 623 return 624 end 625 logical function bsse_create_geom(rtdb) 626 implicit none 627#include "nwc_const.fh" 628#include "errquit.fh" 629#include "inp.fh" 630#include "stdio.fh" 631#include "geom.fh" 632#include "rtdb.fh" 633#include "geomP.fh" 634#include "tcgmsg.fh" 635#include "global.fh" 636#include "mafdecls.fh" 637#include "bsse_common.fh" 638c 639c This logical function generate the geometries of monomers 640c with ghost atoms and store them into DB. 641c 642 integer rtdb ![input] 643 integer geom ![input] 644c 645 character*255 name 646 character*16 tag_old (mx_atm) 647 character*16 tag_new (mx_atm) 648c 649 integer mon_hnd 650 integer k,j,l,m 651c 652 logical is_atm 653 integer atn ![output] 654 double precision mass(mx_atm) 655 character*2 symbol 656 character*16 element 657 658c 659 double precision q_old (mx_atm), q_new (mx_atm) 660 double precision c(3,mx_atm) 661 double precision ctmp(3,mx_atm) 662c 663 logical bsse_rtdb_store 664 external bsse_rtdb_store 665c 666 lenname= inp_strlen(spr_name) 667c 668 if (.not. geom_create(geom,spr_name(1:lenname))) 669 $ call errquit('bsse_create_geom: geom_create failed !', 670 $ 0,GEOM_ERR) 671c 672 if (.not. geom_rtdb_load(rtdb, geom, spr_name)) 673 $ call errquit('bsse_create_geom: no geometry load form rtdb', 674 $ 0,GEOM_ERR) 675c 676 if (.not. geom_cart_get(geom, natoms, tag_old, c, q_old)) 677 $ call errquit('bsse_create_geom: get geom info fail !', 678 $ 0,GEOM_ERR) 679c 680c ------------------------------------------- 681c rename atoms tag to ghost and create geoms 682c ------------------------------------------- 683c 684 do j = 1, nmon 685 do l = 1, natoms 686c 687c initialize all centers as bqX 688c 689 tag_new(l) = 'bq' // tag_old(l) 690 q_new(l) = 0.0d0 691c 692 do k = 1, mon_atm(j) 693c 694c only do bsse with atoms 695c 696 is_atm = geom_tag_to_element(tag_old(l),symbol,element,atn) 697 if ((.not. is_atm) .and. symbol.ne.'bq') 698 $ call errquit('bsse_create_geom: not atom or bq', 699 $ 0,GEOM_ERR) 700c 701c compare and if this is center is an atom then 702c set its tag and charge to the proper values 703c 704 if (l.eq.(mon(j,k))) then 705 tag_new(l) = tag_old(l) 706 q_new(l) = q_old(l) 707 go to 100 708 endif 709 enddo 710 100 continue 711 enddo 712c 713c Sort out the masses of the various centers 714c 715 do l = 1, natoms 716 is_atm = geom_tag_to_element(tag_new(l),symbol,element,atn) 717 if ((.not. is_atm) .and. symbol.ne.'bq') 718 $ call errquit('bsse_create_geom: not atom or bq', 719 $ 0,GEOM_ERR) 720 if (.not.geom_atn_to_default_mass(atn,mass(l))) 721 $ call errquit('bsse_create_geom: default mass failed', 722 $ 911,GEOM_ERR) 723 enddo 724c 725 lenname= inp_strlen(mon_name(j)) 726c 727 if (.not. geom_create(mon_hnd, mon_name(j)(1:lenname)//'g')) 728 $ call errquit('bsse_create_geom: geom_create failed !', 729 $ 0,GEOM_ERR) 730c 731 if (.not. geom_cart_set(mon_hnd, natoms, tag_new, c, q_new)) 732 $ call errquit('bsse_create_geom: geom_cart_set failed', 733 $ 0,GEOM_ERR) 734c 735 if (.not. geom_masses_set(mon_hnd, natoms, mass)) 736 $ call errquit('bsse_create_geom: geom_masses_set failed', 737 $ 0,GEOM_ERR) 738c:debug proposes 739c if (nodeid().eq. 0) then 740c if (.not. geom_print(mon_hnd)) 741c $ call errquit('geom_input: print failed ', 0) 742c endif 743c 744 if (.not.geom_rtdb_store(rtdb, mon_hnd, 745 $ mon_name(j)(1:lenname)//'g')) 746 $ call errquit('bsse_create_geom: geom_rtdb_store failed', 747 $ 60,GEOM_ERR) 748c 749 if (.not. geom_destroy(mon_hnd)) 750 $ call errquit('bsse_create_geom: geom_destroy failed', 751 $ 60,GEOM_ERR) 752c 753c creating and storing monomers without ghost 754c 755 if (.not. geom_create(mon_hnd, mon_name(j)(1:lenname))) 756 $ call errquit('bsse_create_geom: geom_create failed!',0, 757 $ GEOM_ERR) 758c 759 m = 0 760 do l = 1, natoms 761 do k = 1, mon_atm(j) 762 if (l.eq.(mon(j,k))) then 763 m = m + 1 764 tag_new(m) = tag_old(l) 765 q_new (m) = q_old (l) 766 ctmp(1,m) = c(1,l) 767 ctmp(2,m) = c(2,l) 768 ctmp(3,m) = c(3,l) 769 endif 770 enddo 771 enddo 772c 773 if (.not. geom_cart_set(mon_hnd,mon_atm(j),tag_new,ctmp,q_new)) 774 $ call errquit('bsse_create_geom: geom_cart_set failed', 775 $ 0,GEOM_ERR) 776c 777c Work out the masses of the atoms 778c 779 do l=1, mon_atm(j) 780 is_atm = geom_tag_to_element(tag_new(l),symbol,element,atn) 781 if ((.not. is_atm) .and. symbol.ne.'bq') 782 $ call errquit('bsse_create_geom: center is neither atom '// 783 $ 'nor bq',0,GEOM_ERR) 784c.. set default mass 785 if(.not.geom_atn_to_default_mass(atn,mass(l))) 786 & call errquit('bsse_create_geom: default mass failed', 787 & 911,GEOM_ERR) 788 enddo 789c 790 if (.not. geom_masses_set(mon_hnd, mon_atm(j), mass)) 791 $ call errquit('bsse_create_geom: geom_masses_set failed', 792 $ 0,GEOM_ERR) 793c 794 if (nodeid() .eq. 0) then 795 if(.not. geom_print(mon_hnd)) 796 $ call errquit('bsse_create_geom: print failed ', 0, GEOM_ERR) 797 endif 798c 799 if (.not.geom_rtdb_store(rtdb, mon_hnd, mon_name(j)(1:lenname))) 800 $ call errquit('bsse_create_geom: geom_rtdb_store failed', 801 $ 0,GEOM_ERR) 802c 803 if (.not. geom_destroy(mon_hnd)) 804 $ call errquit('bsse_create_geom: geom_destroy failed', 805 $ 0,GEOM_ERR) 806c 807c:debug proposes 808c if (.not. rtdb_print(rtdb, .true.)) call errquit('print failed',0) 809c 810c ====================================================== 811c create basis 812c 813c status=bsse_create_basis(rtdb,geom, tag_new,natoms) 814c======================================================== 815 enddo 816c 817 if(.not. geom_destroy(geom)) 818 $ call errquit('geom_input: geom_destroy failed', 0, GEOM_ERR) 819c 820 bsse_create_geom = .true. 821c 822 return 823 end 824 logical function bsse_energy(rtdb,theory,final_spr_energy) 825 implicit none 826#include "stdio.fh" 827#include "errquit.fh" 828#include "rtdb.fh" 829#include "global.fh" 830#include "mafdecls.fh" 831#include "inp.fh" 832#include "bsse_common.fh" 833c 834 integer rtdb ![input] 835c 836 integer j ! run over monomers 837 integer i 838 integer m_spr 839c tmp 840 character*(*) theory 841 character*255 vec_spr 842 character*255 vec_dbi, vec_dbo 843c 844 character*255 tmp 845c 846 double precision q_spr 847c 848 logical task_energy_doit 849 external task_energy_doit 850c 851 logical bsse_rtdb_load 852 external bsse_rtdb_load 853c 854 logical bsse_create_geom 855 external bsse_create_geom 856c 857c 858 bsse_energy=.false. 859c 860 if(ga_nodeid().eq.0) then 861 write(LuOut,*) 862 write(LuOut,*) 863 call util_print_centered(LuOut, 864 $ 'BSSE Energy Correction', 865 $ 40,.true.) 866 write(LuOut,*) 867 write(LuOut,*) 868 endif 869c 870 if (.not. rtdb_get(rtdb,'bsse:first_j',mt_log,1,first_j)) 871 $ first_j=.true. 872c 873c 874 if(.not.bsse_rtdb_load(rtdb)) 875 $ call errquit('bsse_energy: load data input in db', 911,RTDB_ERR) 876c 877c 878 lentheo = inp_strlen(theory) 879 lenname = inp_strlen(spr_name) 880c 881 if (.not. task_energy_doit(rtdb,theory,spr_energy)) 882 $ call errquit('bsse_energy: no geometry ',0,UNKNOWN_ERR) 883c 884c: take supermolecule total charge 885 if(.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, q_spr)) 886 $ q_spr = 0.0d0 887 888 889c: multiplicity 890c: density methods 891 if ( theory(1:lentheo).eq.'dft' .or. 892 $ theory(1:lentheo).eq.'tddft') then 893 if (.not. rtdb_get(rtdb, 'dft:mult', mt_int, 1, m_spr)) 894 $ call errquit('bsse_energy: rtdb_get of mult failed', 895 $ 0,RTDB_ERR ) 896c: wavefuntion methods 897 898 elseif( theory(1:lentheo).ne.'dft' .and. 899 $ theory(1:lentheo).ne.'tddft') then 900 if (.not. rtdb_get(rtdb,'scf:nopen', MT_INT, 1, m_spr)) 901 $ call errquit('bsse_energy: rtdb_get of nopen failed', 902 $ 0,RTDB_ERR) 903 904 endif 905 906c 907c: name of the original movecs 908 if(theory(1:lentheo).eq.'dft' .or. 909 $ theory(1:lentheo).eq.'tddft') then 910 vec_dbo= 'dft:output vectors' 911 vec_dbi= 'dft:input vectors' 912 elseif (theory(1:lentheo).eq.'mcscf') then 913 vec_dbo= 'mcscf:output vectors' 914 vec_dbi= 'mcscf:input vectors' 915 else 916 vec_dbo = 'scf:output vectors' 917 vec_dbi = 'scf:input vectors' 918 endif 919c 920 if (.not. rtdb_cget(rtdb, vec_dbo, 1, vec_spr)) 921 $ call errquit('bsse_energy: get vectors file failed', 0,RTDB_ERR) 922c 923c: create geom for monomers within supermolecular geom 924 925 926 if(.not.bsse_create_geom(rtdb)) 927 $ call errquit('bsse_energy: bsse_create_geom',911,UNKNOWN_ERR) 928 929 930c 931c:Obtain monomers energies from frozen geometries; 932c:it makes a couple jobs for each monomer (no ghost, ghost) 933c 934 935 j = 1 936 937 do i = 1, nmon*2 938 939 j_mon_name = mon_name(j) 940 lenname = inp_strlen(mon_name(j)) 941c 942 943 if (mod(i,2).eq.0) then 944 945 j_mon_name = j_mon_name(1:lenname)//'g' 946 lenname = lenname + 1 947 948 endif 949c 950 call bsse_param(rtdb, mmon(j), qmon(j), j_mon_name, input(i), 951 $ theory) 952c 953c:evaluate energy 954 if (.not. task_energy_doit(rtdb,theory, mon_energy(i))) 955 $ call 956 $ errquit('bsse_energy: failed calling task_energy',0,UNKNOWN_ERR) 957 958c 959 if (mod(i,2).eq.0) then 960 j = j +1 961 endif 962 963 enddo 964 965 966c----------------------------------------------------------------------- 967c Evaluate BSSE error for supermolecular geometry 968c----------------------------------------------------------------------- 969 970 971 bsse_error = 0.0d0 972 973 i = 1 974 975 do j = 1, nmon 976 977 m_error(j) = mon_energy(i) - mon_energy(i+1) 978 bsse_error = bsse_error + m_error(j) 979 i= i + 2 980 981 enddo 982 983 final_spr_energy = spr_energy + bsse_error 984c 985c:return to original active geom 986 987 988 lenname = inp_strlen(spr_name) 989 990 991 if (.not. rtdb_cput(rtdb,'geometry',1,spr_name(1:lenname))) 992 $ call errquit('bsse_energy: no geometry ',0, RTDB_ERR) 993 994c:return to original output vectors 995c 996 if (.not. rtdb_cput(rtdb, vec_dbi, 1, vec_spr)) 997 & call errquit('bsse_energy: input_vectors',0, RTDB_ERR) 998c 999 1000 if (.not. rtdb_cput(rtdb, vec_dbo, 1,vec_spr)) 1001 & call errquit('bsse_energy: output_vectors',0, RTDB_ERR) 1002 1003 1004c 1005c:return to original charge 1006 1007 1008 if (.not. rtdb_put(rtdb, 'charge', MT_DBL, 1, q_spr)) 1009 $ call errquit 1010 $ ('bsse_energy: failed to write charge to rtdb', 0, RTDB_ERR) 1011 1012 1013c: multiplicity 1014 1015c: density methods 1016 1017 if ( theory(1:lentheo).eq.'dft' .or. 1018 $ theory(1:lentheo).eq.'tddft') then 1019 1020 if (.not. rtdb_put(rtdb, 'dft:mult', mt_int, 1, m_spr)) 1021 $ call errquit('bsse_energy: rtdb_put of mult failed', 1022 $ 0, RTDB_ERR) 1023 1024c: wavefuntion methods 1025 1026 elseif( theory(1:lentheo).ne.'dft' .and. 1027 $ theory(1:lentheo).ne.'tddft') then 1028 1029 if (.not. rtdb_put(rtdb,'scf:nopen', MT_INT, 1, m_spr)) 1030 $ call errquit('bsse_energy: rtdb_put of nopen failed', 1031 $ 0, RTDB_ERR) 1032 1033 1034 endif 1035c 1036c:put into db final energy associated with theory 1037 1038 1039 tmp = theory(1:lentheo)//':energy' 1040c 1041 if (.not. rtdb_put(rtdb,tmp, 1042 $ MT_DBL,1,final_spr_energy)) 1043 $ call 1044 $ errquit('bsse_energy: failed to write charge to rtdb',0,RTDB_ERR) 1045 1046 1047c:debug proposes vama 1048c if (.not. rtdb_print(rtdb, .true.)) call errquit('print failed',0) 1049c 1050c if (nodeid().eq.0) then 1051c do j = 1, nmon 1052c write(LuOut,10) j, m_error(j) 1053c 10 format 1054c $ (/' error for monomer ', i4, ' =', f20.12/) 1055c enddo 1056c 1057 1058 if(ga_nodeid().eq.0) then 1059 1060 write(LuOut,20) bsse_error, spr_energy, final_spr_energy 1061 20 format (/' BSSE error = ', f20.12/ 1062 $ ' Supermolecular energy = ', f20.12/ 1063 $ ' Corrected energy = ', f20.12/) 1064 1065 call util_flush(luout) 1066 1067 endif 1068c 1069 1070 if (.not. rtdb_put(rtdb,'bsse:first_j',mt_log,1,.false.)) 1071 $ call errquit('bsse_energy: rtdb_put failed',0,RTDB_ERR) 1072 1073c 1074 call ga_sync 1075c 1076 bsse_energy = .true. 1077 end 1078 logical function bsse_gradient(rtdb,theory,final_spr_energy, 1079 $ gx_spr) 1080c:debug 1081#include "global.fh" 1082#include "errquit.fh" 1083#include "mafdecls.fh" 1084#include "nwc_const.fh" 1085#include "inp.fh" 1086#include "stdio.fh" 1087#include "bsse_common.fh" 1088#include "rtdb.fh" 1089#include "geom.fh" 1090#include "tcgmsg.fh" 1091c 1092c 1093 integer rtdb ![input] 1094c 1095 integer j,i,l,k,n 1096 integer n_cart 1097 integer n_cartmon 1098 integer geom 1099c 1100 character*255 vec_dbi, vec_dbo 1101 character*255 vec_spr 1102 character*(*) theory 1103 character*16 tag 1104 character*255 tmp 1105c 1106 double precision q_spr 1107 double precision gx_spr(3,*) 1108 double precision gx_mon(3,mx_atm) 1109 double precision crd(3), q 1110c 1111 logical task_gradient_doit 1112 external task_gradient_doit 1113 logical bsse_rtdb_load 1114 external bsse_rtdb_load 1115 logical bsse_create_geom 1116 external bsse_create_geom 1117c 1118c 1119 bsse_gradient=.false. 1120c 1121 if (nodeid().eq.0) then 1122 write(LuOut,*) 1123 write(LuOut,*) 1124 call util_print_centered(LuOut, 1125 $ 'BSSE Correction to Energy Gradient', 1126 $ 40,.true.) 1127 write(LuOut,*) 1128 write(LuOut,*) 1129 endif 1130c 1131 if (.not. rtdb_get(rtdb,'bsse:firt_j',MT_LOG,1,first_j)) 1132 $ first_j=.true. 1133c 1134 if(.not.bsse_rtdb_load(rtdb)) 1135 $ call 1136 $ errquit('bsse_gradient: load data input in db',911,UNKNOWN_ERR) 1137c 1138 1139 lentheo = inp_strlen(theory) 1140 lenname = inp_strlen(spr_name) 1141c 1142 1143 n_cart= 3*natoms 1144c 1145c: get supermolecular gradient and energy 1146 1147 1148 if (.not.task_gradient_doit(rtdb,theory,spr_energy, gx_spr)) 1149 $ call 1150 $ errquit('bsse_gradient: no task gradient do it ',0,UNKNOWN_ERR) 1151 1152 1153c: get supermolecule charge 1154 1155 if (.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, q_spr)) 1156 $ q_spr = 0.0d0 1157 1158 1159c 1160c: multiplicity 1161c 1162c: wavefuntion methods 1163 if( theory(1:lentheo).ne.'dft') then 1164 if (.not. rtdb_get(rtdb,'scf:nopen', MT_INT, 1, m_spr)) 1165 $ call 1166 $ errquit('bsse_gradient: rtdb_get of nopen failed',0,RTDB_ERR) 1167c 1168c: density methods 1169 elseif ( theory(1:lentheo).eq.'dft') then 1170 if (.not. rtdb_get(rtdb, 'dft:mult', mt_int, 1, m_spr)) 1171 $ call 1172 $ errquit('bsse_gradient: rtdb_get failed', 0,RTDB_ERR) 1173 endif 1174c 1175c: name of the original movecs 1176 1177 if(theory(1:lentheo).eq.'dft') then 1178 vec_dbo= 'dft:output vectors' 1179 vec_dbi= 'dft:input vectors' 1180 elseif (theory(1:lentheo).eq.'mcscf') then 1181 vec_dbo= 'mcscf:output vectors' 1182 vec_dbi= 'mcscf:input vectors' 1183 else 1184 vec_dbo = 'scf:output vectors' 1185 vec_dbi = 'scf:input vectors' 1186 endif 1187c 1188 if (.not. rtdb_cget(rtdb, vec_dbo, 1, vec_spr)) 1189 $ call 1190 $ errquit('bsse_gradient: get vectors file failed', 0,RTDB_ERR) 1191c 1192c 1193c: create geom for monomers within supermolecular geom 1194 1195 1196 if(.not.bsse_create_geom(rtdb)) 1197 $ call errquit('bsse_gradient: bsse_create_geom',60,UNKNOWN_ERR) 1198 1199c 1200c: Obtain monomers energies from forzen geometries 1201c: it makes a couple job for each monomer (no ghost, ghost) 1202 j = 1 1203 do i = 1, nmon*2 1204 j_mon_name = mon_name(j) 1205 lenname = inp_strlen(mon_name(j)) 1206 1207c 1208 if (mod(i,2).eq.0) then 1209 j_mon_name = j_mon_name(1:lenname)//'g' 1210 lenname = lenname + 1 1211 endif 1212c 1213 call bsse_param(rtdb, mmon(j), qmon(j), j_mon_name, input(i), 1214 $ theory) 1215c 1216c: evaluate gradient 1217 1218 if (.not.task_gradient_doit(rtdb,theory,mon_energy(i),gx_mon)) 1219 $ call 1220 $ errquit('bsse_gradient: error calling task gradient spr', 1221 $ 0,UNKNOWN_ERR) 1222 1223 if (mod(i,2).eq.0) then 1224 n_cartmon=n_cart 1225 else 1226 n_cartmon=3*mon_atm(j) 1227 endif 1228c 1229c: mix monomers grads with supermolecule grad 1230 if (mod(i,2).eq.0) then 1231 1232 do k = 1, natoms 1233 1234 do n=1 ,3 1235 gx_spr(n,k) = gx_spr(n,k) - gx_mon(n,k) 1236 enddo 1237 1238 enddo 1239c: go for next monomer 1240 j = j +1 1241c 1242 else 1243c 1244 do k = 1 , natoms 1245 do l = 1, mon_atm(j) 1246 if (k.eq.mon(j,l)) then 1247 do n = 1, 3 1248 gx_spr(n,k)=gx_spr(n,k) + gx_mon(n,l) 1249 enddo 1250 endif 1251 enddo 1252 enddo 1253 endif 1254c 1255 enddo 1256 1257 1258c----------------------------------------------------------------------- 1259c Evaluate supermolecular energy free of BSSE 1260c----------------------------------------------------------------------- 1261 1262 bsse_error = 0.0d0 1263 1264 i = 1 1265 1266 do j = 1, nmon 1267 bsse_error = bsse_error + mon_energy(i) - mon_energy(i+1) 1268 i= i + 2 1269 enddo 1270c 1271 1272 final_spr_energy = spr_energy + bsse_error 1273c 1274 1275 lenname = inp_strlen(spr_name) 1276c 1277c 1278 if (.not. rtdb_cput(rtdb,'geometry', 1, spr_name(1:lenname))) 1279 $ call errquit('bsse_gradient: no geometry ',0,RTDB_ERR) 1280 1281 1282c:return to original charge 1283 1284 if (.not. rtdb_put(rtdb, 'charge', MT_DBL, 1, q_spr)) 1285 $ call errquit 1286 $ ('bsse_gradientt: failed to write charge to rtdb', 0,RTDB_ERR) 1287 1288 1289c:return to original output vectors 1290c 1291 if (.not. rtdb_cput(rtdb, vec_dbi, 1, vec_spr)) 1292 & call errquit('bsse_gradient: store input_vectors',60,RTDB_ERR) 1293c 1294 1295 if (.not. rtdb_cput(rtdb, vec_dbo, 1,vec_spr)) 1296 & call errquit('bsse_gradient: store output_vectors',60,RTDB_ERR) 1297 1298 1299 1300c: wavefuntion methods 1301 1302 if( theory(1:lentheo).ne.'dft') then 1303 if (.not. rtdb_put(rtdb,'scf:nopen', MT_INT, 1, m_spr)) 1304 $ call 1305 $ errquit('bsse_gradient: rtdb_put of nopen failed',0,RTDB_ERR) 1306 1307c: density methods 1308 1309 elseif ( theory(1:lentheo).eq.'dft') then 1310 if (.not. rtdb_put(rtdb, 'dft:mult', mt_int, 1, m_spr)) 1311 $ call errquit('bsse_gradient: rtdb_put failed', 0,RTDB_ERR) 1312 endif 1313c 1314c:put into db final gradien associated with theory 1315 1316 1317 tmp=theory(1:lentheo)//':gradient' 1318 if (.not. rtdb_put(rtdb,tmp ,mt_dbl, 1319 $ 3*natoms, gx_spr)) call errquit 1320 $ ('bsse_gradient: could not store gradients',1, RTDB_ERR) 1321c 1322c:put into db final energy associated with theory 1323c 1324 tmp=theory(1:lentheo)//':energy' 1325c 1326 if (.not. rtdb_put(rtdb,tmp, 1327 $ MT_DBL,1,final_spr_energy)) 1328 $ call 1329 $ errquit('bsse_gradient: failed to write charge to rtdb' 1330 $ ,0,RTDB_ERR) 1331 1332c 1333 if (.not. geom_create(geom, spr_name(1:lenname))) call errquit 1334 $ ('bsse_gradient: geom_create final failed !', 0,GEOM_ERR) 1335c 1336c 1337 if (.not.geom_rtdb_load(rtdb, geom, spr_name(1:lenname))) 1338 $ call errquit('bsse_gradient: load geom failed !',0, GEOM_ERR) 1339 1340c:print final energy gradient 1341 1342 if(nodeid().eq.0) then 1343c 1344 write(luout,1000) theory(1:inp_strlen(theory)), 1345 $ 'x','y','z','x','y','z' 1346 1347 do i=1, natoms 1348 if (.not. geom_cent_get(geom, i, tag, crd, q)) call errquit 1349 $ ('bsse_gradient: geometry corrupt?',0, GEOM_ERR) 1350 write(luout,2000) i, tag,(crd(j),j=1,3), 1351 $ (gx_spr(j,i),j=1,3) 1352 enddo 1353 1354 write(luout,*) 1355 write(luout,*) 1356 1357 1000 format(/,/,25X,A,' BSSE ENERGY GRADIENTS',/,/,4X,'atom',15X, 1358 $ 'coordinates', 1359 $ 24X,'gradient',/,6X,2(1X,(3(10X,A1)))) 1360 1361 2000 format(1X,I3,1X,A4,2(1X,3(1X,F10.6))) 1362c 1363c 1364 write(LuOut,20) bsse_error, spr_energy, final_spr_energy 1365 1366 20 format (/' BSSE error = ', f20.12/ 1367 $ ' Supermolecular energy = ', f20.12/ 1368 $ ' Corrected energy = ', f20.12/) 1369 call util_flush(luout) 1370 1371 endif 1372 1373 1374 if(.not. geom_destroy(geom)) 1375 $ call errquit('bsse_gradient: geom_destroy failed', 0,GEOM_ERR) 1376c 1377c 1378 1379 if (.not. rtdb_put(rtdb,'bsse:first_j',mt_log,1,.false.)) 1380 $ call errquit('bsse_gradient: rtdb_put failed',0,RTDB_ERR) 1381c 1382 call ga_sync 1383c 1384 1385 bsse_gradient = .true. 1386 end 1387 logical function bsse_hessian(rtdb) 1388 implicit none 1389#include "util.fh" 1390#include "rtdb.fh" 1391#include "global.fh" 1392#include "mafdecls.fh" 1393#include "inp.fh" 1394#include "bsse_common.fh" 1395#include "stdio.fh" 1396#include "tcgmsg.fh" 1397#include "errquit.fh" 1398c 1399 integer rtdb ![input] 1400c 1401c Call task_hessian_doit for each monomer 1402c and supermolecule and store each one hessian matrix 1403c in different files. As out put write .hess file with 1404c hessian taking account the BSSE 1405c 1406 logical status 1407 integer j, i 1408c 1409 character*255 vec_dbi, vec_dbo 1410 character*255 vec_spr 1411 character*32 theory 1412c 1413 character*16 mult(8) 1414 character*(nw_max_path_len) filehess 1415 character*(nw_max_path_len) filehesstmp 1416 integer m_spr 1417 double precision q_spr 1418 integer lenhess 1419 integer nopen 1420c 1421 logical task_hessian_doit 1422 external task_hessian_doit 1423c 1424 logical bsse_rtdb_load 1425 external bsse_rtdb_load 1426c 1427 logical bsse_create_geom 1428 external bsse_create_geom 1429c 1430c double precision final_spr_energy 1431c 1432 data mult/'singlet','doublet','triplet','quartet', 1433 $ 'quintet','sextet','septet','octet'/ 1434c 1435c call ecce_print_module_entry('task hessian') 1436c 1437 bsse_hessian= .false. 1438c 1439 if (nodeid().eq.0) then 1440 write(LuOut,*) 1441 write(LuOut,*) 1442 call util_print_centered(LuOut, 1443 $ 'BSSE Hessian Correction', 1444 $ 40,.true.) 1445 write(LuOut,*) 1446 write(LuOut,*) 1447 endif 1448c 1449 1450 if (.not. rtdb_get(rtdb,'bsse:firt_j',mt_log,1,first_j)) 1451 $ first_j=.true. 1452c 1453 1454 if(.not.bsse_rtdb_load(rtdb)) 1455 $ call errquit('bsse_hessian: load data input in db', 1456 $ 911,UNKNOWN_ERR ) 1457c 1458c 1459 1460 if (.not. rtdb_cget(rtdb,'task:theory', 1, theory)) 1461 & call errquit('bsse_hessian: get input_vectors',0,RTDB_ERR) 1462 1463 1464 lentheo = inp_strlen(theory) 1465 lenname = inp_strlen(spr_name) 1466c 1467c 1468 1469 if(.not.task_hessian_doit(rtdb)) 1470 $ call errquit('bsse_hessian: error calling hessian doit failed', 1471 $ 0,UNKNOWN_ERR) 1472c 1473 if (.not.rtdb_get 1474 $ (rtdb,'task:energy',mt_dbl,1,spr_energy)) 1475 $ call errquit('bsse_hessian: cannot get energy ',0,RTDB_ERR) 1476 1477 1478c store original sets 1479c: take supermolecule total charge 1480 1481 if(.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, q_spr)) 1482 $ q_spr = 0.0d0 1483 1484c: multiplicity 1485c: wavefuntion methods 1486 1487 if( theory(1:lentheo).ne.'dft') then 1488 1489 if (.not. rtdb_get(rtdb,'scf:nopen', MT_INT, 1, m_spr)) 1490 $ call errquit('bsse_hessian: rtdb_put of nopen failed', 1491 $ 0,RTDB_ERR) 1492 1493c: density methods 1494 1495 elseif ( theory(1:lentheo).eq.'dft') then 1496 1497 if (.not. rtdb_get(rtdb, 'dft:mult', mt_int, 1, m_spr)) 1498 $ call errquit('bsse_hessian: rtdb_put mult failed', 0, RTDB_ERR) 1499 1500 endif 1501c 1502c: name of the original movecs 1503 if(theory(1:lentheo).eq.'dft') then 1504 vec_dbo= 'dft:output vectors' 1505 vec_dbi= 'dft:input vectors' 1506 elseif (theory(1:lentheo).eq.'mcscf') then 1507 vec_dbo= 'mcscf:output vectors' 1508 vec_dbi= 'mcscf:input vectors' 1509 else 1510 vec_dbo = 'scf:output vectors' 1511 vec_dbi = 'scf:input vectors' 1512 endif 1513c 1514 if (.not. rtdb_cget(rtdb, vec_dbo, 1, vec_spr)) 1515 $ call util_file_name('movecs',.false.,.false.,vec_spr) 1516c $ call errquit('bsse_hessian: get vectors file failed', 0) 1517c 1518c: name of hessian file 1519 1520 call util_file_name('hess', .false., .false.,filehess) 1521c:copy to operate the sum at the end| 1522c 1523 if (nodeid().eq.0) then 1524 1525 lenhess= inp_strlen(filehess) 1526 1527 filehesstmp = ' ' 1528 filehesstmp= filehess(1:lenhess)//'.spr' 1529 1530 call util_file_copy(filehess, filehesstmp) 1531 1532 endif 1533 1534c: create geom for monomers within supermolecular geom 1535 1536 if(.not.bsse_create_geom(rtdb)) 1537 $ call errquit('bsse_hessian: bsse_create_geom',911,UNKNOWN_ERR) 1538 1539c 1540c: Obtain monomers energies from forzen geometries; 1541c: it makes a couple jobs for each monomer (no ghost, ghost) 1542c:monomers 1543 1544 j=1 1545 1546 do i=1, nmon*2 1547 1548 j_mon_name = mon_name(j) 1549 lenname = inp_strlen(mon_name(j)) 1550c 1551 if (mod(i,2).eq.0) then 1552 1553 j_mon_name = j_mon_name(1:lenname)//'g' 1554 lenname = lenname + 1 1555 1556 endif 1557c 1558 1559 if (.not. rtdb_cput(rtdb,'geometry', 1, j_mon_name(1:lenname))) 1560 $ call errquit('bsse_hessian: no geometry ',0,RTDB_ERR) 1561c 1562c 1563 1564 call bsse_param(rtdb, mmon(j), qmon(j), j_mon_name, input(i), 1565 $ theory) 1566c 1567c 1568c 1569 if(.not.task_hessian_doit(rtdb)) 1570 $ call errquit('bsse_hessian: error calling hessiandoit', 1571 $ 0,UNKNOWN_ERR) 1572 1573c:get evaluated energy 1574 1575 if (.not.rtdb_get 1576 $ (rtdb,'task:energy', mt_dbl, 1, mon_energy(i))) 1577 $ call errquit('bsse_hessian: get monomer energy',0,RTDB_ERR) 1578 1579c:copy to operate the sum at the end| 1580c 1581 if(nodeid().eq.0) then 1582c 1583 1584 filehesstmp= filehess(1:lenhess)//'.'// 1585 & j_mon_name(1:lenname) 1586 1587 call util_file_copy(filehess, filehesstmp) 1588 1589 endif 1590c 1591 if (mod(i,2).eq.0) then 1592 j = j +1 1593 endif 1594 1595 enddo 1596c 1597 bsse_error = 0.0d0 1598 1599 i = 1 1600 1601 do j = 1, nmon 1602 1603 m_error(j) = mon_energy(i) - mon_energy(i+1) 1604 bsse_error = bsse_error + m_error(j) 1605 i= i + 2 1606 1607 enddo 1608c 1609 final_spr_energy = spr_energy + bsse_error 1610c 1611c sum matrix 1612 1613 if(ga_nodeid().eq.0) then 1614 1615 call hessian_matrix_sum (rtdb) 1616 1617 endif 1618c 1619c: return original geom and parameters 1620c: store supermolecule parameters 1621c 1622 if (.not. rtdb_cput(rtdb, vec_dbi, 1, vec_spr)) 1623 & call errquit('bsse_hessian: put input_vectors',0,RTDB_ERR) 1624c 1625 1626 if (.not. rtdb_cput(rtdb, vec_dbo, 1,vec_spr)) 1627 & call errquit('bsse_hessian: output_vectors',0,RTDB_ERR) 1628 1629 1630 1631c:return to original active geom 1632 1633c 1634 lenname = inp_strlen(spr_name) 1635 1636 if (.not. rtdb_cput(rtdb,'geometry', 1, spr_name(1:lenname))) 1637 $ call errquit ('bsse_hessian: no geometry ',0,RTDB_ERR) 1638 1639c: store supermolecule total charge 1640 1641 if (.not. rtdb_get(rtdb, 'charge', MT_DBL, 1, q_spr)) 1642 $ q_spr = 0.0d0 1643 1644c: store multiplicity 1645c: wavefuntion methods 1646 1647 if( theory(1:lentheo).ne.'dft') then 1648 1649 if (.not. rtdb_put(rtdb,'scf:nopen', MT_INT, 1, m_spr)) 1650 $ call errquit('bsse_hessian: rtdb_put of nopen failed' 1651 $ , 0, RTDB_ERR) 1652 1653c: density methods 1654 1655 elseif ( theory(1:lentheo).eq.'dft') then 1656 1657 if (.not. rtdb_put(rtdb, 'dft:mult', mt_int, 1, m_spr)) 1658 $ call errquit('bsse_gradient:cant put m_spr rtdb',0, RTDB_ERR) 1659 1660 endif 1661 1662c 1663 1664 if(ga_nodeid().eq.0) then 1665 1666 write(LuOut,20) bsse_error, spr_energy, final_spr_energy 1667 20 format (/' BSSE error = ', f20.12/ 1668 $ ' Supermolecular energy = ', f20.12/ 1669 $ ' Corrected energy = ', f20.12/) 1670 1671 call util_flush(luout) 1672 1673 endif 1674 1675 1676 if (.not. rtdb_put(rtdb,'bsse:first_j',mt_log,1,.false.)) 1677 $ call errquit('bsse_hessian: rtdb_put failed',0, RTDB_ERR) 1678 1679c 1680 1681 1682 bsse_hessian = .true. 1683 1684c 1685 call ga_sync 1686c 1687 end 1688 subroutine hessian_matrix_sum(rtdb) 1689c 1690c This function makes the apropiate sum to mix the hessian of monomers matrixes 1691C (with ghost and without ghost) to supermolecular hessian matrix 1692c 1693 implicit none 1694#include "global.fh" 1695#include "mafdecls.fh" 1696#include "rtdb.fh" 1697#include "util.fh" 1698#include "inp.fh" 1699#include "stdio.fh" 1700#include "bsse_common.fh" 1701#include "errquit.fh" 1702c 1703 integer rtdb 1704 integer j,i 1705 integer iii 1706 integer k 1707 integer nat2 1708 integer lenhess 1709 double precision dval 1710 1711 integer nhesst 1712 character*(nw_max_path_len) filehess 1713c 1714 integer k_exy_mont 1715 integer k_exy_mon 1716 integer k_exy_spr 1717 integer k_exycp 1718 integer k_exybt 1719 integer k_exybs 1720c 1721 integer l_exycp 1722 integer l_exybt 1723 integer l_exy_mon 1724 integer l_exy_mont 1725 integer l_exybs 1726 integer n3xyz 1727c 1728 n3xyz=3*natoms 1729 nhesst=n3xyz*(n3xyz+1)/2 1730 nat2=n3xyz*n3xyz 1731c 1732 if (.not.ma_push_get(mt_dbl,nat2,'bsse exybs ', 1733 * l_exybs,k_exybs)) 1734 * call errquit('bsse_hessian_exy: cannot allocate',650, MA_ERR) 1735 call dfill(nat2,0.0d00,dbl_mb(k_exybs),1) 1736c 1737 if (.not.ma_push_get(mt_dbl,nhesst,'bsse exybt ', 1738 * l_exybt,k_exybt)) 1739 * call errquit('bsse_hessian_exy: cannot allocate',650, MA_ERR) 1740 call dfill(nhesst,0.0d00,dbl_mb(k_exybt),1) 1741c 1742 if (.not.ma_push_get(mt_dbl,nat2,'bsse exycp ', 1743 * l_exycp,k_exycp)) 1744 * call errquit('bsse_hessian_exy: cannot allocate',650, MA_ERR) 1745 call dfill(nat2,0.0d00,dbl_mb(k_exycp),1) 1746c 1747 if (.not.ma_push_get(mt_dbl,nat2,'bsse:hessian:exy mon', 1748 * l_exy_mon,k_exy_mon)) 1749 * call errquit('bsse_hessian_exy: cannot allocate',650, MA_ERR) 1750c 1751 if (.not.ma_push_get(MT_DBL,nhesst,'bsse:hessian:exyt mon', 1752 * l_exy_mont,k_exy_mont)) 1753 * call errquit('bsse_hessian_exy: cannot allocate',650, MA_ERR) 1754c 1755c:read triangle matrix super and do square 1756 1757 1758 call util_file_name('hess', .false., .false.,filehess) 1759c 1760 1761 lenhess=inp_strlen(filehess) 1762 filehess=filehess(1:lenhess)//'.spr' 1763c 1764 1765 open(unit=69,file=filehess,form='formatted',status='old', 1766 & err=99900,access='sequential') 1767 do iii = 0,(nhesst-1) 1768 read(69,*,err=99901,end=99902) dval 1769 dbl_mb(k_exybt+iii) = dval 1770 enddo 1771 close(unit=69,status='keep') 1772c 1773 1774 call trin2squa(dbl_mb(k_exybs),dbl_mb(k_exybt),n3xyz) 1775c 1776 1777 j=1 1778 1779 do i=1, nmon*2 1780 1781 j_mon_name = mon_name(j) 1782 lenname = inp_strlen(mon_name(j)) 1783 1784 if (mod(i,2).eq.0) then 1785 1786 j_mon_name = j_mon_name(1:lenname)//'g' 1787 1788 lenname = lenname + 1 1789 filehess=filehess(1:lenhess)//'.'//j_mon_name 1790 1791 n3xyz=3*natoms 1792 1793 nhesst = n3xyz*(n3xyz+1)/2 1794c:read triangle matrix monomer and do square 1795 open(unit=69,file=filehess,form='formatted',status='old', 1796 & err=99900,access='sequential') 1797 1798 do iii = 0,(nhesst - 1) 1799 read(69,*,err=99901,end=99902) dval 1800 1801 dbl_mb(k_exy_mont+iii) = dval 1802 1803 enddo 1804 1805 close(unit=69,status='keep') 1806 1807 call trin2squa(dbl_mb(k_exy_mon),dbl_mb(k_exy_mont),n3xyz) 1808c:substract entire ghost matrix case 1809 do k = 0,(nat2- 1) 1810 1811 dval=dbl_mb(k_exybs+k)-dbl_mb(k_exy_mon+k) 1812 dbl_mb(k_exybs+k)=dval 1813 1814 enddo 1815c:visit other monomer 1816 1817 j = j +1 1818 1819 else 1820 1821 n3xyz=3*mon_atm(j) 1822 1823 nhesst = n3xyz*(n3xyz+1)/2 1824c:read triangle matrix monomer 1825 1826 filehess=filehess(1:lenhess)//'.'//j_mon_name 1827c 1828 1829 open(unit=69,file=filehess,form='formatted',status='old', 1830 1831 & err=99900,access='sequential') 1832 1833 do iii = 0,(nhesst - 1) 1834 1835 read(69,*,err=99901,end=99902) dval 1836 dbl_mb(k_exy_mont+iii) = dval 1837 1838 enddo 1839 1840 close(unit=69,status='keep') 1841 1842c: do square 1843 call trin2squa(dbl_mb(k_exy_mon),dbl_mb(k_exy_mont),n3xyz) 1844 1845c:sum of no-ghost monomer case 1846 call dcopy(nat2,dbl_mb(k_exybs),1,dbl_mb(k_exycp),1) 1847 1848c 1849 call sum_small_matrix(dbl_mb(k_exybs),dbl_mb(k_exycp), 1850 $ dbl_mb(k_exy_mon),3*natoms, n3xyz, j) 1851 endif 1852 enddo 1853c 1854c n3xyz=3*natoms 1855c 1856c:write the final hessian 1857c 1858 call util_file_name('hess', .false., .false.,filehess) 1859c 1860 call stpr_wrt_fd_from_sq(dbl_mb(k_exybs),n3xyz,filehess) 1861 1862c:clean memory 1863 1864 if (.not.ma_chop_stack(l_exy_mont)) 1865 * call errquit('hess_read: cannot deallocate Exyt',555, MA_ERR) 1866 1867 if (.not.ma_chop_stack(l_exy_mon)) 1868 * call errquit('hess_read: cannot deallocate Exyt',555, MA_ERR) 1869 1870 if (.not.ma_chop_stack(l_exycp)) 1871 * call errquit('hess_read: cannot deallocate l_Exycp',555, MA_ERR) 1872 1873 if (.not.ma_chop_stack(l_exybt)) 1874 * call errquit('hess_read: cannot deallocate l_Exybt',555, MA_ERR) 1875 1876 if (.not.ma_chop_stack(l_exybs)) 1877 * call errquit('hess_read: cannot deallocate l_Exybs',555, MA_ERR) 1878c 1879 return 188099900 continue 1881 write(luout,*)'hess_file => ',filehess 1882 call errquit('error opening file: "filehess"', 1883 1 DISK_ERR,911) 188499901 continue 1885 write(luout,*)'hess_file => ',filehess 1886 call errquit('error reading file: "filehess"', 1887 1 DISK_ERR,911) 188899902 continue 1889 write(luout,*)'hess_file => ',filehess 1890 call errquit('unexpected EOF reading file: "filehess"', 1891 1 DISK_ERR,911) 1892 end 1893c 1894 subroutine sum_small_matrix(d,a,b,dim_a,dim_b,j) 1895c 1896c this routine sums a monomer matrix in fractions 3x3 of each atom 1897c to the supermolecular matrix 1898c vama 1899#include "bsse_common.fh" 1900 integer i,k,m,n,l,g,r,c 1901 integer dim_a ! dimension matrix super 1902 integer dim_b ! dimension matrix monomer 1903 integer j ! monomer j, used to take the list of atoms 1904 double precision d(dim_a,dim_a) ! matrix super resul 1905 double precision a(dim_a,dim_a) ! matrix super 1906 double precision b(dim_b,dim_b) ! matrix monomer 1907 double precision dval 1908c 1909 do i = 1 , natoms 1910 do k = 1 , natoms 1911c 1912 do m= 1, mon_atm(j) 1913 do n= 1, mon_atm(j) 1914c 1915 l= mon(j,n) 1916 g= mon(j,m) 1917c 1918 if(l.eq.k.and.g.eq.i) then 1919c 1920 do r=0, 2 1921 do c=0,2 1922 dval=a((i-1)*3+1+r,(k-1)*3+1 +c)+ 1923 $ b((m-1)*3+1+r,(n-1)*3+1+c) 1924 1925 d((i-1)*3+1+r,(k-1)*3+1+c)=dval 1926 d((k-1)*3+1+c,(i-1)*3+1+r)=dval 1927 enddo 1928 enddo 1929 endif 1930 enddo 1931 enddo 1932 enddo 1933 enddo 1934 end 1935c 1936 subroutine trin2squa(a,b,dime) 1937c 1938c convert a triangular matrix to a simetric square matrix 1939c 1940 integer dime 1941 integer j, l, i 1942 double precision a(dime,dime) !squ 1943 double precision b((dime*(dime+1))/2) !tri 1944 double precision dval 1945c 1946 i = 1 1947 do j=1, dime 1948 do l= 1, j 1949 dval=b(i) 1950 a(j,l)=dval 1951 a(l,j)=dval 1952 i=i+1 1953 enddo 1954 enddo 1955 end 1956 subroutine cfill(n,val,a,ia) 1957 implicit none 1958 integer n, ia 1959 character*(*) val, a(*) 1960 integer i 1961c 1962c initialise characters array 1963c 1964 if (ia.eq.1) then 1965 do 10 i = 1, n 1966 a(i) = val 1967 10 continue 1968 else 1969 do 20 i = 1,(n-1)*ia+1,ia 1970 a(i) = val 1971 20 continue 1972 endif 1973c 1974 end 1975 1976c $Id$ 1977