1 logical function drv_lst(rtdb) 2* 3* $Id$ 4* 5 implicit none 6#include "mafdecls.fh" 7#include "rtdb.fh" 8#include "global.fh" 9#include "stdio.fh" 10#include "util.fh" 11 integer rtdb 12 logical hnd_lstx 13 external hnd_lstx 14 integer ir, iw 15 logical some, dbug 16 common/hnd_iofile/ir,iw 17c 18 call util_print_push 19 call util_print_rtdb_load(rtdb, 'driver') 20 call ecce_print_module_entry('driver') 21 if (util_print('brdcst', print_never)) call setdbg(1) 22c 23 dbug=.false. 24 some=.false. 25 some=some.or.dbug 26 ir=LuIn 27 iw=LuOut 28 if(some) then 29 write(iw,*) ' drv_lst: calling hnd_lstx . ' 30 endif 31c 32c ----- lst path ----- 33c 34 drv_lst = hnd_lstx(rtdb) 35c 36 if (drv_lst) then 37 call ecce_print_module_exit('driver', 'ok') 38 else 39 call ecce_print_module_exit('driver', 'failed') 40 endif 41c 42 if (util_print('brdcst', print_never)) call setdbg(0) 43 call util_print_pop 44c 45 end 46 logical function hnd_lstx(rtdb) 47 implicit double precision (a-h,o-z) 48#include "errquit.fh" 49#include "rtdb.fh" 50#include "msgtypesf.h" 51#include "global.fh" 52#include "mafdecls.fh" 53#include "stdio.fh" 54#include "geom.fh" 55#include "util.fh" 56 logical geom_zmt_get_nzvar 57 logical geom_zmt_get_nizmat 58 external geom_zmt_get_nizmat 59 logical geom_zmt_get_nzfrz 60 external geom_zmt_get_nzfrz 61 logical geom_zmt_get_izfrz 62 external geom_zmt_get_izfrz 63 logical geom_lst_get_coord 64 external geom_lst_get_coord 65 integer rtdb 66 integer geom 67 logical do_gradient 68 logical some 69 logical out 70 logical dbug 71 logical status 72 logical done 73 logical rstart 74 logical zcoord 75 character*16 atmlab 76 parameter (mxatom=500) 77 parameter (mxcart=3*mxatom) 78 parameter (mxzmat=1500) 79 parameter (mxcoor=1500) 80c mxcoor=max(mxcart,mxzmat) 81 common/hnd_iofile/ir,iw 82 common/hnd_optmiz/x0(mxcoor),x(mxcoor),dx(mxcoor), 83 1 g0(mxcoor),g(mxcoor),ds(mxcoor), 84 2 func,func0,gmax,gmax0,curv,alpha,gnorm 85 common/hnd_optfrz/nzfrz,izfrz(mxcoor),iatfrz(mxatom) 86 common/hnd_optvar/zcoord,ncoord,mcoord 87 common/hnd_optfun/e,eg(mxcart) 88 common/hnd_lstrun/nptlst,iptlst 89 common/hnd_lsttim/energy_time,gradient_time 90 common/hnd_zmtpar/nzmat,nzvar,nvar 91 common/hnd_molxyz/c(mxcart),zan(mxatom),nat 92 common/hnd_mollab/atmlab(mxatom) 93 data zero /0.0d+00/ 94 data pt5 /0.5d+00/ 95 data one /1.0d+00/ 96c 97 hnd_lstx = .false. 98c 99 do_gradient = .false. 100c 101 dbug=.false. 102 out =.false. 103 out =out.or.dbug 104 some=.false. 105 some=some.or.out 106 if (ga_nodeid().eq.0.or.some) then 107 write(iw,9999) 108 endif 109c 110c ----- get going ... ----- 111c 112 if (.not. geom_create(geom, 'geometry')) 113 & call errquit('hnd_lst: geom_create?', 911, GEOM_ERR) 114 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 115 & call errquit('hnd_lst: no geometry ', 911, RTDB_ERR) 116 if (.not. geom_cart_get(geom, ncent, atmlab, c, zan)) 117 & call errquit('hnd_lst: geom_get ', 911, GEOM_ERR) 118 nat = ncent 119 if(.not.geom_zmt_get_nizmat(geom,nzmat)) 120 & call errquit('geom_input: geom_zmt_get_nizmat failed',0, 121 & GEOM_ERR) 122 if(.not.geom_zmt_get_nzvar(geom,nzvar)) 123 & call errquit('geom_input: geom_zmt_get_nzvar failed',0, 124 & GEOM_ERR) 125c 126c ----- scan -lst- coordinates file ----- 127c 128 nptlst=1 129 10 continue 130 status=geom_lst_get_coord(x(1),x(1+nat),x(1+nat), 131 & nat,nptlst) 132 if(status) then 133 nptlst=nptlst+1 134 go to 10 135 else 136 nptlst=nptlst-1 137 endif 138c 139 if (ga_nodeid().eq.0.or.dbug) then 140 write(iw,9996) nptlst 141 endif 142 if(nptlst.eq.0) then 143 write(iw,*) 'no -lst- points found. return and continue ..' 144 return 145 endif 146c 147 nzfrz=0 148 do i=1,mxcoor 149 izfrz(i)=0 150 enddo 151 do i=1,mxatom 152 iatfrz(i)=0 153 enddo 154c 155c ----- frozen coordinates ? if so .... ----- 156c 157 nzfrz=0 158 if (geom_zmt_get_nzfrz(geom,nzfrz)) then 159 if(.not.geom_zmt_get_izfrz(geom,izfrz,nzfrz)) 160 $ call errquit('geom_input: geom_zmt_get_izfrz failed',0, 161 & GEOM_ERR) 162 endif 163 if (ga_nodeid().eq.0.or.dbug) then 164 if(nzfrz.gt.0) then 165 write(iw,9993) nzfrz 166 write(iw,9992) (izfrz(i),i=1,nzfrz) 167 else 168 write(iw,9991) 169 endif 170 endif 171c 172c ----- get ready to start now ----- 173c 174c if internal coordinates ... 175c nzvar = # of (redundant) internal coordinates 176c ncoord = # of independent internal coordinates 177c mcoord = max( 3*nat , nzvar ) 178c if cartesian coordinates ... 179c ncoord = 3*nat 180c mcoord = 3*nat 181c 182 rstart =.false. 183c 184 zcoord=.false. 185 zcoord=nzmat.gt.0 186 if(dbug) then 187 write(iw,*) 'nzmat,zcoord = ',nzmat,zcoord 188 write(iw,*) 'nzfrz = ',nzfrz 189 if(nzfrz.gt.0) then 190 write(iw,*) 'izfrz = ',(izfrz(i),i=1,nzfrz) 191 endif 192 endif 193 if(zcoord) then 194 call geom_bandbi(geom) 195 ncart =3*nat 196 ncoord=nvar 197 mcoord=max(ncart,nzvar) 198 else 199 nzmat =0 200 nzvar =0 201 nvar =3*nat 202 ncart =3*nat 203 ncoord=ncart 204 mcoord=ncart 205 endif 206 if(dbug) then 207 write(iw,9997) nzvar,nvar,ncart,ncoord,mcoord 208 endif 209c 210c ----- set up ---- 211c 212 e =zero 213 func =zero 214 func0 =zero 215 gmax =zero 216 gmax0 =zero 217 curv =zero 218 do i=1,mxcoor 219 dx(i)=zero 220 ds(i)=zero 221 eg(i)=zero 222 g0(i)=zero 223 x0(i)=zero 224 x(i)=zero 225 g(i)=zero 226 enddo 227c 228c ----- first point ----- 229c 230 iptlst=1 231 do i=1,ncart 232 x(i)=c(i) 233 enddo 234c 235c ----- energy ----- 236c 237 time_start=util_wallsec() 238 call hnd_lst_energy(rtdb,geom) 239 time_end =util_wallsec() 240 energy_time=time_end-time_start 241 if(out) then 242 write(iw,*) 'energy_time = ',energy_time 243 endif 244c 245 func=e 246 if(do_gradient) then 247c 248c ----- gradient ----- 249c 250 time_start=util_wallsec() 251 call hnd_opt_gradient(rtdb,geom) 252 time_end =util_wallsec() 253 gradient_time=time_end-time_start 254 if(out) then 255 write(iw,*) 'gradient_time = ',gradient_time 256 endif 257c 258 do i=1,ncart 259 g(i)=eg(i) 260 enddo 261 if(zcoord) then 262 call hnd_opt_tfgx(rtdb,geom) 263 endif 264 endif 265c 266 call hnd_lst_print(rtdb) 267 do i=1,ncart 268 x0(i)=x(i) 269 enddo 270 do i=1,ncoord 271 g0(i)=g(i) 272 enddo 273c 274 if (ga_nodeid().eq.0.and.dbug) then 275 if (.not. geom_print(geom)) call errquit 276 $ ('hnd_lst_drv: geom_print?',0, GEOM_ERR) 277 if (util_print('bonds',print_default)) then 278 if (.not.geom_print_distances(geom)) call errquit( 279 & 'hnd_lst_drv: geom_print_distances failed',911, 280 & GEOM_ERR) 281 endif 282 if (util_print('angles',print_default)) then 283 if (.not.geom_print_angles(geom)) call errquit( 284 & 'hnd_lst_drv: geom_print_angles failed',911, GEOM_ERR) 285 endif 286 endif 287 call ga_sync() 288c 289c ----- lst path ----- 290c 291 100 continue 292 iptlst=iptlst+1 293c 294 call hnd_lst_path(rtdb,geom,rstart) 295c 296 if(zcoord) then 297 call hnd_opt_tfgx(rtdb,geom) 298 endif 299 call hnd_lst_print(rtdb) 300c 301 if (ga_nodeid().eq.0.and.dbug) then 302 if (.not. geom_print(geom)) call errquit 303 $ ('hnd_lst_drv: geom_print?',0, GEOM_ERR) 304 if (util_print('bonds',print_default)) then 305 if (.not.geom_print_distances(geom)) call errquit( 306 & 'hnd_lst_drv: geom_print_distances failed',911, GEOM_ERR) 307 endif 308 if (util_print('angles',print_default)) then 309 if (.not.geom_print_angles(geom)) call errquit( 310 & 'hnd_lst_drv: geom_print_angles failed',911, 311 & GEOM_ERR) 312 endif 313 endif 314 call ga_sync() 315c 316c ----- go back to next -lst- point ----- 317c 318 done=iptlst.eq.nptlst 319 if(.not.done) then 320 go to 100 321 else 322 if (ga_nodeid().eq.0.or.dbug) then 323 write(iw,9998) 324 endif 325 call ga_sync() 326 if (.not.geom_destroy(geom)) 327 & call errquit('hnd_lst: geom_destroy?',911, GEOM_ERR) 328 hnd_lstx=.true. 329 endif 330c 331 9999 format(/,10x,8(1h-), 332 1 /,10x,'LST Path', 333 2 /,10x,8(1h-)) 334 9998 format(/,10x,18(1h-), 335 1 /,10x,'LST Path completed', 336 2 /,10x,18(1h-)) 337 9997 format(' in hnd_lstx, nzvar,nvar,ncart,ncoord,mcoord = ',5i5) 338 9996 format(' number of points on -lst- path (nptlst) = ',i4) 339 9993 format(' number of frozen internal coordinates = ',i4) 340 9992 format(' frozen internal coordinates = ',12i4) 341 9991 format(' no frozen internal coordinates.') 342 end 343 SUBROUTINE HND_LST_PRINT(RTDB) 344 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 345#include "mafdecls.fh" 346#include "rtdb.fh" 347#include "global.fh" 348 INTEGER RTDB 349 PARAMETER (MXATOM=500) 350 PARAMETER (MXCART=3*MXATOM) 351 PARAMETER (MXZMAT=1500) 352 PARAMETER (MXCOOR=1500) 353 LOGICAL DBUG 354 CHARACTER*16 ATMNAM 355 COMMON/HND_IOFILE/IR,IW 356 COMMON/HND_MOLLAB/ATMNAM(MXATOM) 357 COMMON/HND_MOLNUC/NUC(MXATOM) 358 COMMON/HND_MOLXYZ/C(3,MXATOM),ZAN(MXATOM),NAT 359 COMMON/HND_OPTMIZ/X0(MXCOOR),X(MXCOOR),DX(MXCOOR), 360 1 G0(MXCOOR),G(MXCOOR),DS(MXCOOR), 361 2 FUNC,FUNC0,GMAX,GMAX0,CURV,ALPHA,GNORM 362 COMMON/HND_LSTRUN/NPTLST,IPTLST 363 COMMON/HND_OPTFRZ/NZFRZ,IZFRZ(MXCOOR),IATFRZ(MXATOM) 364 COMMON/HND_OPTFUN/E,EG(MXCART) 365C 366 DBUG=.FALSE. 367C 368C ----- PRINT OPTIMIZATION SUMMARY ----- 369C 370 IF( GA_NODEID().EQ.0.OR.DBUG) THEN 371 WRITE(IW,9999) NPTLST,IPTLST,FUNC 372 ENDIF 373 RETURN 374 9999 FORMAT(1H1,/,1X, 375 1 'nptlst iptlst func ',1X,I5,I8,F17.8) 376 END 377 SUBROUTINE HND_LST_PATH(RTDB,GEOM,RSTART) 378 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 379#include "mafdecls.fh" 380#include "tcgmsg.fh" 381#include "msgtypesf.h" 382#include "global.fh" 383#include "geom.fh" 384#include "rtdb.fh" 385C 386C ----- ONE DIMENSIONAL SEARCH. ----- 387C 388 PARAMETER (MXATOM=500) 389 PARAMETER (MXCART=3*MXATOM) 390 PARAMETER (MXZMAT=1500) 391 PARAMETER (MXCOOR=1500) 392C 393 INTEGER RTDB 394 INTEGER GEOM 395 LOGICAL STATUS 396 LOGICAL UTIL_TEST_TIME_REMAINING 397 EXTERNAL UTIL_TEST_TIME_REMAINING 398 LOGICAL GEOM_LST_GET_COORD 399 EXTERNAL GEOM_LST_GET_COORD 400 CHARACTER*16 TAGS_NW 401 DIMENSION COORDS_NW(MXCART) 402 DIMENSION CHARGE_NW(MXATOM) 403 DIMENSION TAGS_NW(MXATOM) 404C 405 LOGICAL DO_GRADIENT 406 LOGICAL RSTART 407 LOGICAL DBUG 408 LOGICAL OUT 409 COMMON/HND_IOFILE/IR,IW 410 COMMON/HND_MOLXYZ/C(MXCART),ZAN(MXATOM),NAT 411 COMMON/HND_OPTMIZ/X0(MXCOOR),X(MXCOOR),DX(MXCOOR), 412 1 G0(MXCOOR),G(MXCOOR),DS(MXCOOR), 413 2 FUNC,FUNC0,GMAX,GMAX0,DELSQ,ALPHA,GNORM 414 COMMON/HND_LSTRUN/NPTLST,IPTLST 415 COMMON/HND_OPTFRZ/NZFRZ,IZFRZ(MXCOOR),IATFRZ(MXATOM) 416 COMMON/HND_OPTTIM/ENERGY_TIME,GRADIENT_TIME 417 COMMON/HND_OPTFUN/E,EG(MXCART) 418 COMMON/HND_ZMTPAR/NZMAT,NZVAR,NVAR 419 DATA ZERO,TWO,THREE /0.0D+00,2.0D+00,3.0D+00/ 420 DATA ONEPT5 /1.5D+00/ 421C 422 DO_GRADIENT=.FALSE. 423C 424 DBUG=.FALSE. 425 OUT =.FALSE. 426 OUT =OUT.OR.DBUG 427 IF(OUT) THEN 428 WRITE(IW,9994) 429 ENDIF 430C 431 NCART =3*NAT 432C 433C ----- NORMAL SETUP ----- 434C 435 FUNC0 = FUNC 436 DO I = 1,NCART 437 X0(I) = X(I) 438 G0(I) =EG(I) 439 G (I) =EG(I) 440 X (I) =ZERO 441 ENDDO 442C 443C ----- NEW -LST- POINT ----- 444C 445 STATUS=GEOM_LST_GET_COORD(X(1),X(1+NAT),X(1+NAT*2), 446 1 NAT,IPTLST) 447 DO IAT = 1,NAT 448 C(1+3*(IAT-1)) = X(IAT) 449 C(2+3*(IAT-1)) = X(IAT+NAT) 450 C(3+3*(IAT-1)) = X(IAT+NAT*2) 451 ENDDO 452 DO I = 1,NCART 453 X(I)=ZERO 454 X(I)=C(I) 455 ENDDO 456C 457C ----- WRITE TO -NWCHEM- ----- 458C 459 STATUS=GEOM_CART_GET(GEOM,NAT_NW,TAGS_NW,COORDS_NW, 460 1 CHARGE_NW) 461C 462 IF(DBUG) THEN 463 WRITE(IW,9997) 464 WRITE(IW,9998) ( X0(I),I=1,NCART) 465 WRITE(IW,9996) 466 WRITE(IW,9998) ( X(I),I=1,NCART) 467 WRITE(IW,9999) 468 WRITE(IW,9998) (COORDS_NW(I),I=1,NCART) 469 ENDIF 470 DO I = 1,NCART 471 COORDS_NW(I)=X(I) 472 ENDDO 473 STATUS=GEOM_CART_SET(GEOM,NAT_NW,TAGS_NW,COORDS_NW, 474 1 CHARGE_NW) 475 STATUS=GEOM_RTDB_STORE(RTDB,GEOM,'geometry') 476C 477 IF(OUT) THEN 478 WRITE(IW,9999) 479 WRITE(IW,9998) (COORDS_NW(I),I=1,NCART) 480 ENDIF 481C 482C 483C ----- CALL FUNCTION EVALUATION ----- 484C 485 STATUS=UTIL_TEST_TIME_REMAINING(RTDB,INT(ENERGY_TIME*ONEPT5)) 486 IF(OUT) THEN 487 WRITE(IW,*) 488 1 'ENOUGH TIME REMAINING FOR ENERGY = ? ',STATUS 489 ENDIF 490 IF(.NOT.STATUS) THEN 491 WRITE(IW,*) 492 1 'NOT ENOUGH TIME REMAINING, SHUTTING DOWN ... ' 493 RETURN 494 ENDIF 495C 496 CALL HND_LST_ENERGY(RTDB,GEOM) 497 FUNC=E 498C 499C ----- CALCULATE THE GRADIENT FOR THE FINAL POINT ----- 500C 501 900 CONTINUE 502 IF(DO_GRADIENT) THEN 503 STATUS=UTIL_TEST_TIME_REMAINING(RTDB,INT(GRADIENT_TIME*ONEPT5)) 504 IF(OUT) THEN 505 WRITE(IW,*) 506 1 'ENOUGH TIME REMAINING FOR GRADIENT = ? ',STATUS 507 ENDIF 508 IF(.NOT.STATUS) THEN 509 WRITE(IW,*) 510 1 'NOT ENOUGH TIME REMAINING, SHUTTING DOWN ... ' 511 RETURN 512 ENDIF 513C 514 CALL HND_OPT_GRADIENT(RTDB,GEOM) 515 DO I=1,NCART 516 G(I) =EG(I) 517 ENDDO 518 ENDIF 519C 520C ----- RETURN AFTER SUCCESSFUL POINT ----- 521C 522 1000 CONTINUE 523 IF(DBUG) THEN 524 WRITE(IW,9993) 525 ENDIF 526 RETURN 527 9999 FORMAT(' IN SEARCH_LIN, COORDS = ') 528 9998 FORMAT(F12.7) 529 9997 FORMAT(' IN SEARCH_LIN, X0 = ') 530 9996 FORMAT(' IN SEARCH_LIN, DX = ') 531 9995 FORMAT(' IN SEARCH_LIN, X = ') 532 9994 FORMAT(' STARTING SEARCH_LIN ') 533 9993 FORMAT(' ENDING SEARCH_LIN ') 534 END 535 subroutine hnd_lst_energy(rtdb,geom) 536 implicit double precision (a-h,o-z) 537#include "errquit.fh" 538#include "mafdecls.fh" 539#include "rtdb.fh" 540#include "geom.fh" 541#include "global.fh" 542 integer rtdb 543 integer geom 544 logical task_energy 545 external task_energy 546 logical status 547 logical dbug 548 parameter (mxatom=500) 549 parameter (mxcart=3*mxatom) 550 parameter (mxzmat=1500) 551 parameter (mxcoor=1500) 552 common/hnd_iofile/ir,iw 553 common/hnd_lstrun/nptlst,iptlst 554 common/hnd_optfun/e,eg(mxcart) 555 common/hnd_molxyz/c(3,mxatom),zan(mxatom),nat 556 common/hnd_zmtpar/nzmat,nzvar,nvar 557c 558 logical geom_print_zmatrix 559c 560 dbug=.false. 561c 562 if( ga_nodeid().eq.0.or.dbug) then 563 write(iw,9999) nptlst,iptlst 564 if(dbug) then 565 write(iw,9996) 566 do iat=1,nat 567 write(iw,9995) iat,(c(i,iat),i=1,3) 568 enddo 569 write(iw,9997) 570 endif 571 if(.not.geom_print(geom)) 572 1 call errquit('hnd_lst_energy: print error',911, GEOM_ERR) 573 endif 574c 575 if(nzmat.gt.0) then 576 call geom_bandbi(geom) 577 if (.not. geom_print_zmatrix(geom,0d0,.false.)) 578 $ call errquit('hnd_lst: print zmat failed',0, GEOM_ERR) 579 write(iw,9997) 580 endif 581c 582 if (task_energy(rtdb)) then 583 status=rtdb_get(rtdb,'task:energy',MT_DBL,1,e) 584 if(dbug) then 585 write(iw,9998) e 586 endif 587 else 588 call errquit('lst path: energy failed', 0, GEOM_ERR) 589 endif 590c 591 return 592 9999 format(1h1,' nptlst',i3,//,' iptlst',I3) 593 9998 format(' in lst_energy, e = ',f15.10) 594 9997 format(/) 595 9996 format(21x,21(1h-),/, 596 1 21x,'cartesian coordinates',/, 597 2 21x,21(1h-)) 598 9995 format(9x,i5,3f15.8) 599 end 600