1* 2* $Id$ 3* 4* this file contains auxilary routines to stpr_gen functions 5* 6* current routines: 7* stpr_fd_upd_dipole ! computes finite difference dipole moment 8* stpr_fd_upd_hess ! computes finite difference either central or forward 9* stpr_wrt_fd_from_sq ! writes hessian to file 10* stpr_check_genat_restart ! check for restart "is info available?" 11* stpr_get_genat_restart ! get restart info 12* stpr_put_genat_restart ! put restart info out to restart file 13* stpr_gen_hess_foldave ! averages off diaginal contributions 14* stpr_gen_hess_fold ! sums off diaginal contributions (partial computations) 15* stpr_gen_set_diag ! sets inactive atom diagonal contribs to large value 16* 17c 18C> \ingroup stpr_priv 19C> @{ 20C> 21 subroutine stpr_fd_upd_dipole(ddipole,mdipole,pdipole, 22 & s_delta,delta,nat,iatom,ixyz,q1) 23 implicit none 24c 25c::passed 26 integer nat 27 integer iatom 28 integer ixyz 29 double precision s_delta 30 double precision delta 31 double precision mdipole(3) 32 double precision pdipole(3) 33 double precision ddipole(3,3,nat) 34 double precision q1 35c:local 36 integer moment 37 double precision rdelta 38 double precision value 39c 40 rdelta = 1.0d00/(s_delta*delta) 41 do moment = 1,3 42 value = rdelta*(pdipole(moment)-mdipole(moment)) * q1 43 ddipole(moment,ixyz,iatom) = value 44 enddo 45c 46 end 47C> 48C> \brief Update the Hessian from new gradient data 49C> 50C> Apply the finite difference formula to the gradient data provided 51C> to calculate a contribution to and update the Hessian. 52C> 53C> Depending on the values of the parameters `gradm`, `gradp` and 54C> `s_delta` the subroutine will evaluate either the forward difference 55C> equation 56C> \f{eqnarray*}{ 57C> f'(x) &=& \frac{f(x+h)-f(x)}{h} 58C> \f} 59C> or the central difference formula 60C> \f{eqnarray*}{ 61C> f'(x) &=& \frac{f(x+h)-f(x-h)}{2h} 62C> \f} 63C> 64 subroutine stpr_fd_upd_hess(rtdb, 65 & hess,gradm,gradp,s_delta,delta,nat, 66 & iatom_t,ixyz_t) 67 implicit none 68#include "errquit.fh" 69#include "stdio.fh" 70#include "geom.fh" 71#include "sym.fh" 72c::passed 73 integer rtdb !< [Input] The RTDB handle 74 integer nat !< [Input] The number of atoms 75 integer iatom_t !< [Input] The displaced atom 76 integer ixyz_t !< [Input] The displaced coordinate 77 double precision hess(3,nat,3,nat) !< [In/Output] The Hessian 78 double precision gradm(3,nat) !< [Input] The gradient at 79 !< negative (or no) displacement 80 double precision gradp(3,nat) !< [Input] The gradient at 81 !< positive displacement 82 double precision delta !< [Input] The step size 83 double precision s_delta !< [Input] The scale factor 84 !< for total displacement 85 !< * 1 - for forward difference 86 !< * 2 - for central difference 87c::local 88 integer geom 89 integer iatom,ixyz 90 double precision rdelta, value 91 double precision q2 92c 93 if (.not.geom_create(geom,'reference')) call errquit 94 $ ('stpr_fd_upd_hess:geom_create failed?',1, GEOM_ERR) 95 if (.not.geom_rtdb_load(rtdb,geom,'reference')) call errquit 96 $ ('stpr_fd_upd_hess:geom_rtdb_load failed?',2, RTDB_ERR) 97c 98c finite difference [g(x+delta) - g(x-delta)]/(s_delta*delta) (s_delta = 2.0) 99c central difference [g(x+delta) - g(x)]/(s_delta*delta) (s_delta = 1.0) 100c 101c 102 rdelta = 1.0d00/(s_delta*delta) 103 do 00100 iatom = 1,iatom_t 104 if (sym_atom_pair(geom,iatom_t,iatom,q2)) then 105** write(6,*) ' iatom_t iatom q2 ', iatom_t, iatom, q2 106 if (iatom.ne.iatom_t) q2 = q2 + q2 107 do 00200 ixyz = 1,3 108 value = rdelta*(gradp(ixyz,iatom)-gradm(ixyz,iatom)) 109 value = q2*value 110 hess(ixyz_t,iatom_t,ixyz,iatom) = value 111** hess(ixyz,iatom,ixyz_t,iatom_t) = value 11200200 continue 113 endif 11400100 continue 115c 116 if (.not.geom_destroy(geom)) 117 $ call errquit 118 $ ('stpr_fd_upd_hess: geom_destroy failed?',33, GEOM_ERR) 119 end 120 subroutine stpr_wrt_fd_dipole(ddipole,nat,filename) 121 implicit none 122#include "stdio.fh" 123#include "inp.fh" 124#include "util.fh" 125 integer nat 126 double precision ddipole(3,3,nat) 127 character*(*) filename 128c 129 integer lu 130 integer print_level 131 integer atom, xyz, moment 132 logical does_it_exist 133c 134 call util_print_get_level(print_level) 135 lu = 67 136 does_it_exist = .false. 137 inquire(file=filename,exist=does_it_exist) 138 if ((does_it_exist).and.(print_level.gt.print_none)) 139 & write(luout,*) 140 & 'stpr_wrt_fd_dipole: overwrite of existing file', 141 & filename(1:inp_strlen(filename)) 142 open(unit=lu,file=filename, 143 & form='formatted', 144 & access='sequential', 145 & status='unknown') 146c 147 do atom = 1,nat 148 do xyz = 1,3 149 do moment = 1,3 150 write(lu,10000)ddipole(moment,xyz,atom) 151 enddo 152 enddo 153 enddo 154c 15510000 format(1x,1pd20.10) 156c 157 close(unit=lu,status='keep') 158c 159 end 160 subroutine stpr_wrt_fd_from_sq(hess,rank_hess,filename) 161 implicit none 162c 163#include "stdio.fh" 164#include "inp.fh" 165#include "util.fh" 166c 167 integer rank_hess 168 double precision hess(rank_hess,rank_hess) 169 character*(*) filename 170c 171 logical does_it_exist 172c 173 integer i, j, lu 174 integer print_level 175c 176 call util_print_get_level(print_level) 177 lu = 66 178 does_it_exist = .false. 179 inquire(file=filename,exist=does_it_exist) 180 if ((does_it_exist).and.(print_level.gt.print_none)) 181 & write(luout,*) 182 & ' stpr_wrt_fd_from_sq: overwrite of existing file:', 183 & filename(1:inp_strlen(filename)) 184 open(unit=lu,file=filename, 185 & form='formatted', 186 & access='sequential', 187 & status='unknown') 188c 189 do 00100 i = 1,rank_hess 190 do 00200 j = 1,i 191 write(lu,10000)hess(i,j) 19200200 continue 19300100 continue 194c 19510000 format(1x,1pd20.10) 196c 197 close(unit=lu,status='keep') 198c 199 end 200C> 201C> \brief Check and retrieve restart data 202C> 203C> Check the runtime database for the state of the calculation 204C> and report back whether this calculation needs to be restarted 205C> or not. The function also provides the rank of the initial 206C> atom and atomic coordinate to start from. 207C> 208C> The initial atom and atomic coordinate are read from the 209C> restart file. For this purpose the function may open the restart 210C> file if available. But the restart file must be closed before 211C> the function returns. 212C> 213C> \return Returns .TRUE. if the calculation needs to be restarted, 214C> and .FALSE. otherwise. 215C> 216 logical function stpr_check_genat_restart( 217 & rtdb, iatom_start,ixyz_start) 218 implicit none 219#include "errquit.fh" 220#include "global.fh" 221#include "mafdecls.fh" 222#include "msgtypesf.h" 223#include "msgids.fh" 224#include "stdio.fh" 225#include "cstprfiles.fh" 226 integer rtdb !< [Input] The RTDB handle 227 integer iatom_start !< [Output] The atom to start from 228 integer ixyz_start !< [Output] The atomic coordinate to start 229 !< from 230c 231 integer int_restart 232 integer rank, ijunk1, ijunk2 233 integer mitob1 234 logical does_it_exist 235c 236 logical ostart, orestart, ocontinue 237c 238c assume not a restart 239 iatom_start = 1 240 ixyz_start = 1 241 int_restart = 0 242 call util_get_rtdb_state(rtdb,ostart,ocontinue,orestart) 243 if (ostart) then 244 if (ga_nodeid().eq.0) 245 & call util_file_unlink(FILEATR) 246 else if (ocontinue.or.orestart) then 247 if (ga_nodeid().eq.0) then 248 does_it_exist = .false. 249 inquire(file=FILEATR,exist=does_it_exist) 250 if (does_it_exist) then 251 open(unit=69,file=FILEATR, 252 & form='unformatted', 253 & access='sequential', 254 & status='old') 255 read(69)iatom_start,ixyz_start,rank, ijunk1, ijunk2 256 close(unit=69,status='keep') 257 int_restart = 1 258 else 259 write(luout,*)'*** Warning continue called for but no ***' 260 write(luout,*)'*** fd restart file for nuclear hessian ***' 261 write(luout,*)'*** starting from scratch so to speak ***' 262 endif 263 endif 264 else 265 call errquit 266 & ('stpr_check_genat_restart: error with rtdb state',911, 267 & RTDB_ERR) 268 endif 269c 270 mitob1=MA_sizeof(MT_INT,1,MT_BYTE) 271 call ga_brdcst(Msg_gen_at_iatom +MSGINT,iatom_start,mitob1,0) 272 call ga_brdcst(Msg_gen_at_ixyz +MSGINT,ixyz_start, mitob1,0) 273 call ga_brdcst(Msg_gen_at_restart+MSGINT,int_restart,mitob1,0) 274c 275 if (int_restart.eq.1) then 276 stpr_check_genat_restart = .true. 277 else if (int_restart.eq.0) then 278 stpr_check_genat_restart = .false. 279 else 280 write(luout,*)' invalid int_restart value ', ga_nodeid() 281 call errquit(' stpr_check_genat_restart: fatal error ', 282 & int_restart, INPUT_ERR) 283 endif 284 end 285C> 286C> \brief Load the restart information 287C> 288C> Load the data needed for restarting a Hessian calculation. 289C> 290C> Obviously the restart data is read from the restart file. Hence, 291C> this subroutine must open the restart file, but the restart file 292C> must be closed before the routine returns. 293C> 294 subroutine stpr_get_genat_restart(rank_in,hess,grad0,get_grad0, 295 & dipole_okay,ddipole) 296 implicit none 297#include "errquit.fh" 298#include "global.fh" 299#include "stdio.fh" 300#include "util.fh" 301#include "cstprfiles.fh" 302 integer rank_in !< [Input] The dimension of the Hessian 303 double precision hess(rank_in,rank_in) !< [Output] The Hessian 304 double precision grad0(rank_in) !< [Output] The gradient 305 double precision ddipole(3*rank_in) !< [Output] The dipole 306 !< derivative 307 logical get_grad0 !< [Input] Should the gradient be loaded 308 !< * `.TRUE.` - load the gradient 309 !< * `.FALSE.` - do not load the gradient 310 logical dipole_okay !< [Output] Is the dipole derivative returned 311 !< * `.TRUE.` - The dipole derivative is 312 !< returned 313 !< * `.FALSE.` - The dipole derivative is 314 !< not returned 315c 316 logical does_it_exist 317 integer ijunk1, ijunk2, rank, iflag_grad0 318 integer dipole_there 319c 320 if (ga_nodeid().ne.0) then 321 write(luout,*)' non-master node called me ',ga_nodeid() 322 call errquit('stpr_get_genat_restart: error ',911, INPUT_ERR) 323 endif 324c 325 inquire(file=FILEATR,exist=does_it_exist) 326 if (does_it_exist) then 327 open(unit=69,file=FILEATR, 328 & form='unformatted', 329 & access='sequential', 330 & status='old') 331 read(69)ijunk1,ijunk2,rank,iflag_grad0,dipole_there 332 if (dipole_there.eq.1) then 333 dipole_okay = .true. 334 else 335 dipole_okay = .false. 336 endif 337 if (rank.ne.rank_in) then 338 write(luout,*)'rank not the same as rank_in ' 339 write(luout,*)' rank :',rank 340 write(luout,*)' rank_in :',rank_in 341 close(unit=69,status='keep') 342 call errquit('stpr_get_genat_restart: error ',911, INPUT_ERR) 343 endif 344 if (get_grad0.and.iflag_grad0.ne.1) then 345 write(luout,*)' grad 0 not written but requested ' 346 call errquit(' stpr_get_genat_restart: error',911, INPUT_ERR) 347 endif 348 if ((.not.get_grad0).and.iflag_grad0.eq.1) then 349 write(luout,*)' grad 0 written but not requested ' 350 call errquit(' stpr_get_genat_restart: error',911, INPUT_ERR) 351 endif 352 if (get_grad0) read(69) grad0 353 read(69) hess 354 if (dipole_okay) read(69) ddipole 355 if (util_print('debug_stepper_restart',print_debug) 356 & .or. 357 & util_print('debug_stepper',print_debug)) then 358 write(6,*)'hessian read from restart file ' 359 call output(hess,1,rank,1,rank,rank,rank,1) 360 call stpr_print_ddipole(ddipole, 361 & 'dipole derivative read from restart file', 362 & (rank/3), 363 & 1.0d-07) 364 endif 365 close(unit=69,status='keep') 366 else 367 write(6,*)' no finite difference hessian restart ', 368 & 'information read ' 369 endif 370 end 371C> 372C> \brief Store the restart data 373C> 374C> Write the restart data to a file. The restart file is opened and 375C> overwritten with new restart data. The restart file must be closed 376C> before the subroutine returns. 377C> 378 subroutine stpr_put_genat_restart(rank,hess,grad0, 379 & iatom_in,ixyz_in,nat,put_grad0, 380 & dipole_okay,ddipole) 381 implicit none 382#include "errquit.fh" 383#include "global.fh" 384#include "stdio.fh" 385#include "util.fh" 386#include "cstprfiles.fh" 387 integer rank !< [Input] The dimension of the Hessian 388 integer iatom_in !< [Input] The current atom 389 integer ixyz_in !< [Input] The current atomic coordinate 390 integer nat !< [Input] The number of atoms 391 double precision hess(rank,rank) !< [Input] The Hessian 392 double precision grad0(rank) !< [Input] The gradient 393 double precision ddipole(3*rank) !< [Input] The dipole derivative 394 logical put_grad0 !< [Input] Should the gradient be stored? 395 logical dipole_okay !< [Input] Should the dipole derivative be 396 !< stored? 397 logical lopen !< Is the file open? 398c 399 integer iatom_start, ixyz_start 400 integer iflag_grad0 401 integer dipole_there 402c 403 if (ga_nodeid().ne.0) then 404 write(luout,*)' non-master node called me ',ga_nodeid() 405 call errquit('stpr_put_genat_restart: error ',911, INPUT_ERR) 406 endif 407c 408 if(iatom_in.eq.nat.and.ixyz_in.eq.3) then 409 call util_file_unlink(FILEATR) 410 return 411 endif 412 iatom_start = iatom_in 413 ixyz_start = ixyz_in + 1 414 if(ixyz_in.eq.3) then 415 iatom_start = iatom_start + 1 416 ixyz_start = 1 417 endif 418c 419 dipole_there = 0 420 if (dipole_okay) dipole_there = 1 421c 422 if (put_grad0) then 423 iflag_grad0 = 1 424 else 425 iflag_grad0 = 0 426 endif 427c If unit 69 is not closed then the OPEN statement will be ignored 428c in which case the restart data is appended to the existing file 429c rather than overwriting it. 430c Therefore we force the file to be closed. If it is not then there 431c is a bug somewhere else. 432 inquire(unit=69,opened=lopen) 433 if (lopen) call errquit("stpr_put_genat_restart: unit 69 should " 434 & //"be closed at this point",0,UERR) 435 open(unit=69,file=FILEATR, 436 & form='unformatted', 437 & access='sequential', 438 & status='unknown') 439 write(69)iatom_start,ixyz_start,rank,iflag_grad0,dipole_there 440 if (put_grad0) write(69)grad0 441 write(69)hess 442 if (dipole_okay)write(69) ddipole 443 close(unit=69,status='keep') 444 if (util_print('debug_stepper_restart',print_debug) 445 & .or. 446 & util_print('debug_stepper',print_debug)) then 447 write(6,*)'hessian put to restart file ' 448 call output(hess,1,rank,1,rank,rank,rank,1) 449 call stpr_print_ddipole(ddipole, 450 & 'dipole derivative put to restart file', 451 & (rank/3), 452 & 1.0d-07) 453 endif 454 end 455 subroutine stpr_gen_hess_foldave(hess,rank_hess) 456*! averages off diaginal contributions 457 implicit none 458 integer rank_hess 459 double precision hess(rank_hess,rank_hess) 460* 461 integer i,j 462 double precision dbl_tmp 463* 464 do i = 1,rank_hess 465 do j = 1,(i-1) 466 dbl_tmp = hess(i,j) + hess(j,i) 467 dbl_tmp = dbl_tmp/2.0d00 468 hess(i,j) = dbl_tmp 469 hess(j,i) = dbl_tmp 470 enddo 471 enddo 472 end 473 subroutine stpr_gen_hess_fold(hess,rank_hess) 474*! sums off diaginal contributions assuming a partial computation 475 implicit none 476#include "util.fh" 477 integer rank_hess 478 double precision hess(rank_hess,rank_hess) 479* 480 integer i,j 481 double precision dbl_tmp 482 integer icount 483 double precision dbl_diff, max_dbl_diff 484 logical o_debug 485c 486 o_debug = util_print('debug_stepper_restart',print_debug) 487 o_debug = o_debug .or. 488 & util_print('debug_stepper',print_debug) 489 if (o_debug) then 490 write(6,*)' hessian before fold operaton' 491 call output(hess,1,rank_hess,1,rank_hess, 492 & rank_hess,rank_hess,1) 493 icount = 0 494 max_dbl_diff = -1.0d00 495 endif 496* 497 do i = 1,rank_hess 498 do j = 1,(i-1) 499 if (o_debug) then 500 dbl_diff = abs(hess(i,j)) - abs(hess(j,i)) 501 max_dbl_diff = max(max_dbl_diff,dbl_diff) 502 icount = icount + 1 503 write(6,12345)icount,dbl_diff,max_dbl_diff 504 endif 505 dbl_tmp = hess(i,j) + hess(j,i) 506 hess(i,j) = dbl_tmp 507 hess(j,i) = dbl_tmp 508 enddo 509 enddo 510 if (o_debug) then 511 write(6,*)' hessian after fold operaton' 512 call output(hess,1,rank_hess,1,rank_hess, 513 & rank_hess,rank_hess,1) 514 endif 51512345 format('<',i2,'> <diff=',f14.8,'> <diff_max=',f14.8) 516 end 517 subroutine stpr_gen_set_diag(hess,rank_hess) 518*! sets diag to default value for stiff frequency analysis 519*! e.g., active atom computation 520 implicit none 521#include "util.fh" 522 integer rank_hess 523 double precision hess(rank_hess,rank_hess) 524* 525 integer i 526 double precision dbl_tmp 527 logical o_debug 528 o_debug = util_print('debug_stepper_restart',print_debug) 529 o_debug = o_debug .or. 530 & util_print('debug_stepper',print_debug) 531 if (o_debug) then 532 write(6,*)' hessian before diag set operaton' 533 call output(hess,1,rank_hess,1,rank_hess, 534 & rank_hess,rank_hess,1) 535 endif 536* 537 dbl_tmp = 1.0d00 538 do i = 1,rank_hess 539 if (hess(i,i).eq.0.0d00) hess(i,i) = dbl_tmp 540 enddo 541 if (o_debug) then 542 write(6,*)' hessian after diag set operaton' 543 call output(hess,1,rank_hess,1,rank_hess, 544 & rank_hess,rank_hess,1) 545 endif 546 end 547 subroutine stpr_print_ddipole(ddipole,msg,nat,thresh) 548 implicit none 549#include "errquit.fh" 550#include "stdio.fh" 551#include "inp.fh" 552 integer nat 553 double precision ddipole(3,3,nat) 554 double precision thresh 555 character*(*) msg 556 character*1 cname(3) 557 character*1 kapname 558c 559 double precision val 560 integer moment 561 integer xyz, atom 562c 563 cname(1) = 'x' 564 cname(2) = 'y' 565 cname(3) = 'z' 566c 567 write(luout,'(/a/)') msg(1:inp_strlen(msg)) 568 write(luout,*)' ' 569 do moment = 1,3 570 kapname = cname(moment) 571 call inp_ucase(kapname) 572 write(luout,*)' ' 573 write(luout,10001) kapname 574 do atom = 1,nat 575 do xyz = 1,3 576 if (moment.ge.1.and.moment.le.3) then 577 val = ddipole(moment,xyz,atom) 578* 579* from A Physicists Desk Reference, The Second Edition of Physics Vade Mecum 580* Herbert L. Anderson, Editor in Chief 581* Copyright (C) 1989 American Institute of Physics 582* 335 East 45th Street, New York, NY 10017 583* 584*1 debye = 10**(-18) esu cm * [1 e/4.8032068 x 10**(-10) esu]*[1 m /100cm]*[a0/5.29177249 x 10**(-11) m] 585*1 debye = (1.0/4.8032068/5.29177249) * 10**(-18 + 10 - 2 + 11) e a0 586*1 debye = (1.0/4.8032068/5.29177249) * 10**(1) e a0 587*1 e a0 = (4.8032068*5.29177249) * 10**(-1) debye 588*1 e a0 = 25.417477608 * 10**(-1) debye 589*1 e a0 = 2.5417477608 debye 590* 591*use 1 e a0 = 2.541 7478 debye 592* 593 val = val*2.5417478d00 594 val = val/0.529177249d00 ! bohr->angstrom matches current geom data 595 if (abs(val).gt.thresh) then 596 write(luout,10002) 597 & cname(moment),atom,cname(xyz), 598 & ddipole(moment,xyz,atom),val 599 endif 600 else 601 write(luout,10000)moment 602 call errquit('stpr_print_ddipole: fatal error',911, 603 & UNKNOWN_ERR) 604 endif 605 enddo 606 enddo 607 enddo 608 write(luout,*)' ' 609 write(luout,*)' ' 61010000 format('invalid moment value',i10) 61110001 format(1x,a1,1x, 612 & 'vector of derivative dipole (au) [debye/angstrom]') 61310002 format(1x,'d_dipole_',a1,'/<atom=',i4,',',a1,'> = ', 614 & f10.4,5x,'[',f10.4,']') 615 end 616C> 617C> @} 618