1C> \ingroup task 2C> @{ 3 logical function task_mepgs(rtdb) 4 implicit none 5#include "errquit.fh" 6#include "global.fh" 7#include "mafdecls.fh" 8#include "stdio.fh" 9#include "util.fh" 10#include "nwc_const.fh" 11#include "cgsopt.fh" 12#include "cmepgs.fh" 13#include "rtdb.fh" 14#include "geom.fh" 15c 16 integer rtdb 17c 18 integer geom 19 logical ignore 20 logical mepgs_freq 21 logical mepgs_opt 22 logical firstpass 23 double precision energyts 24 logical status 25 double precision start ! Tracks time used in last step 26 logical task_gradient, task_energy 27 external task_gradient, task_energy 28c 29c Disable printing to ecce of movecs after this point 30c 31 call movecs_ecce_print_off() 32 firstpass = .true. 33c 34 1000 continue 35c 36 stotal = 0d0 37c 38c Read input, load /cmepgs/, get geometry 39c 40 call mepgs_initialize(rtdb, geom, firstpass) 41c 42c **** Obtain initial energy **** 43c 44 if (firstpass) then 45 if (.not. task_energy(rtdb)) 46 $ call errquit('mepgs: task_energy failed',0, GEOM_ERR) 47 if (.not. rtdb_get(rtdb,'task:energy', mt_dbl, 1, energy)) 48 $ call errquit('mepgs: could not get energy',0, RTDB_ERR) 49 energyts = energy 50 end if 51c 52c **** Print trayectory file **** 53c 54 energy = energyts 55 call mepgs_path(geom, .true., .false.) 56ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 57c Displace TS along selected mode ! has to be negative c 58c Assumes the user is sure on the TS nature and that a c 59c preferently a frequency analysis has already been performed c 60ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 61c 62c **** Deallocate geom **** 63c 64 if (.not.geom_destroy(geom)) 65 & call errquit('hnd_opt: geom_destroy?', 911, GEOM_ERR) 66c 67c **** Move away from TS ***** 68c 69 if (firstpass) then 70 if (.not. mepgs_freq(rtdb)) 71 $ call errquit('mepgs: mepgs_freq failed',0, GEOM_ERR) 72 end if 73c 74c *** Reallocate geom info **** 75c 76 if (.not. geom_create(geom, 'geometry')) 77 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 78c 79c **** Select side to traverse **** 80c 81 if (forward) then 82 if (.not. geom_rtdb_load(rtdb, geom, 'ircforward')) 83 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 84 if (ga_nodeid().eq.0) write(6,5000) 855000 format(/,10x,24('-'),/,10x,'Forward IRC optimization', 86 $ /,10x,24('-')) 87 else if (backward) then 88 if (.not. geom_rtdb_load(rtdb, geom, 'ircbackward')) 89 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 90 if (ga_nodeid().eq.0) write(6,5100) 915100 format(/,10x,25('-'),/,10x,'Backward IRC optimization', 92 $ /,10x,25('-')) 93 end if 94c 95c **** Store selected side **** 96c 97 if (.not. geom_rtdb_store(rtdb, geom, 'geometry')) 98 $ call errquit('gsopt_energy_step: grs?',geom, RTDB_ERR) 99 if (.not. geom_ncent(geom,nat)) 100 $ call errquit('hnd_opt: natoms?',nat, GEOM_ERR) 101 call grad_active_atoms(rtdb, nat, oactive, nactive) 102 if (.not. geom_systype_get(geom, isystype)) 103 $ call errquit('mepgs: systype?',0, GEOM_ERR) 104c 105c **** Energy and Gradient **** 106c 107 call mepgs_gra(rtdb, geom) 108 if (.not. rtdb_get(rtdb,'task:energy', mt_dbl, 1, energy)) 109 $ call errquit('mepgs: could not get energy',0, RTDB_ERR) 110c 111c **** Construct projector **** 112c 113 call gsopt_cart_pmat(rtdb, geom) 114 call dcopy(ncart, gx, 1, g, 1) ! g() set to gx() 115c 116c **** Print trayectory file **** 117c 118 call mepgs_path(geom, .false., .false.) 119c 120c **** Check energy decrease agreement **** 121c 122 if (ga_nodeid().eq.0) write(6,6000) evib 1236000 format(/,10x,'Expected decrease in energy', 10x, f12.6) 124 if (ga_nodeid().eq.0) write(6,6100) energyts - energy 1256100 format(/,10x,'Obtained decrease in energy', 10x, f12.6) 126 if ((energyts - energy) .lt. 0.0d0) then 127 if (ga_nodeid().eq.0) write(6,6200) 1286200 format(/,1x,25('-')/,1x,'The energy has increased'/,1x,25('-') ) 129 ircdone = .false. 130 goto 2000 131 end if 132c 133c **** Obtain displacement step **** 134c 135 call gsopt_compute_actual_step(geom) 136c 137c **** Read Initial hessian for IRC and OPT loops **** 138c **** to be able to handle separate updates **** 139c 140 call mepgs_hss_init(rtdb,geom) 141c 142c **** Update the IRC hessian **** 143c 144CCCCCCCCCCCCCCCCCCCCCC 145CCCC MASS CCCC 146CCCCCCCCCCCCCCCCCCCCCC 147 if (mswg) then 148 call mwcoord( ds, nvar, .true.) 149 call mwgrad( g, nvar, .true.) 150 call mwgrad(oldgra, nvar, .true.) 151 end if 152CCCCCCCCCCCCCCCCCCCCCC 153CCCC MASS CCCC 154CCCCCCCCCCCCCCCCCCCCCC 155 call mepgs_hessian_update() 156CCCCCCCCCCCCCCCCCCCCCC 157CCCC MASS CCCC 158CCCCCCCCCCCCCCCCCCCCCC 159 if (mswg) then 160 call mwcoord( ds, nvar, .false.) 161 call mwgrad( g, nvar, .false.) 162 call mwgrad(oldgra, nvar, .false.) 163 end if 164CCCCCCCCCCCCCCCCCCCCCC 165CCCC MASS CCCC 166CCCCCCCCCCCCCCCCCCCCCC 167c 168c **** initialization **** 169c 170 flip = .false. 171 ircdone = .false. 172c 173c **** Deallocate geom **** 174c 175 if (.not.geom_destroy(geom)) 176 & call errquit('hnd_opt: geom_destroy?', 911, GEOM_ERR) 177c 178c **** Gonzalez & Schlegel Iterative loop **** 179c 180 if (.not. mepgs_opt(rtdb)) 181 $ call errquit('mepgs: could not optimize',0, RTDB_ERR) 182c 183c **** Obtain second side **** 184c 185 if (ircboth) then 186 forward = .false. 187 ircboth = .false. 188 backward = .true. 189 firstpass = .false. 190 goto 1000 191 end if 192c 193 2000 continue 194c 195 task_mepgs=ircdone 196c 197 call ga_sync() 198c 199 end 200C> @} 201CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc 202CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc 203 subroutine mepgs_initialize(rtdb, geom, start) 204 implicit none 205#include "errquit.fh" 206#include "nwc_const.fh" 207#include "cgsopt.fh" 208#include "cmepgs.fh" 209#include "geom.fh" 210#include "rtdb.fh" 211#include "util.fh" 212#include "global.fh" 213#include "mafdecls.fh" 214#include "inp.fh" 215 integer rtdb 216 integer geom ! [output] 217 logical start ! [output] 218c 219c This routine initializes the common /cmepgs/ and 220c also creates and returns the geometry handle 221c 222 integer i, j, num, ma_type, nactive_atoms, l_actlist 223 logical ignore 224 character*80 title 225 character*8 source, test 226 character*32 theory 227 logical gsopt_geom_cart_coords_get 228c 229 call util_print_push 230 call util_print_rtdb_load(rtdb, 'mepgs') 231 call ecce_print_module_entry('mepgs') 232 oprint = util_print('information', print_low) 233 $ .and. (ga_nodeid() .eq. 0) 234 odebug = util_print('debug', print_debug) 235 $ .and. (ga_nodeid() .eq. 0) 236c 237 if (rtdb_cget(rtdb,'title',1,title)) then 238 if (oprint) then 239 write(6,*) 240 write(6,*) 241 call util_print_centered(6, title, 40, .false.) 242 write(6,*) 243 write(6,*) 244 endif 245 endif 246c 247c ----- parameters for optimization mepgs ----- 248c 249 if (.not. rtdb_get(rtdb,'mepgs:evib',mt_dbl,1,evib)) 250 $ evib = 0.0001d0 251 if (.not. rtdb_get(rtdb,'mepgs:stride',mt_dbl,1,stride)) 252 $ stride = 0.1d0 253 if (.not. rtdb_cget(rtdb,'mepgs:xyz',1,xyz)) 254 $ xyz = ' ' 255 if (.not. rtdb_get(rtdb,'mepgs:inhess',mt_int,1,inhess)) 256 $ inhess=0 257 if (.not. rtdb_get(rtdb,'mepgs:nircopt',mt_int,1,nircopt)) 258 $ nircopt=250 259 if (.not. rtdb_get(rtdb,'ircgs:mswg',mt_log,1,mswg)) 260 $ mswg = .false. 261c 262c **** Select side to traverse **** 263c 264 if (start) then 265 ircboth = .true. 266 forward = .true. 267 if (rtdb_get(rtdb,'mepgs:backward',mt_log,1,backward)) then 268 backward = .true. 269 forward = .false. 270 ircboth = .false. 271 end if 272 if (rtdb_get(rtdb,'mepgs:forward',mt_log,1,forward)) then 273 backward = .false. 274 forward = .true. 275 ircboth = .false. 276 end if 277c 278c Save a copy of the TS geometry 279c 280 if (ga_nodeid() .eq. 0) then 281 ignore = rtdb_parallel(.false.) 282 if (.not. geom_create(geom, 'geometry')) 283 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 284 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 285 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 286 if (.not. geom_rtdb_store(rtdb, geom, 'tsreference')) 287 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 288 if (.not. geom_destroy(geom)) 289 $ call errquit('mepgs: geom_destroy?',0, GEOM_ERR) 290 ignore = rtdb_parallel(.true.) 291 endif 292 end if 293 call ga_sync() 294c 295c Load the geometry info 296c 297 if (.not. geom_create(geom, 'geometry')) 298 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 299 if (start) then 300 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 301 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 302 else if (.not. start) then 303 if (.not. geom_rtdb_load(rtdb, geom, 'tsreference')) 304 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 305 end if 306CJMC mass 307 if (start .and. mswg) then 308 if (.not. geom_masses_get(geom, nat, atmass)) 309 & call errquit('ircgs: geom_masses_get failed',911, GEOM_ERR) 310 end if 311CJMC mass 312c 313 if (.not. geom_ncent(geom,nat)) 314 $ call errquit('hnd_opt: natoms?',nat, GEOM_ERR) 315 call grad_active_atoms(rtdb, nat, oactive, nactive) 316 if (.not. geom_systype_get(geom, isystype)) 317 $ call errquit('mepgs: systype?',0, GEOM_ERR) 318c 319 if (oprint) then 320 if (ga_nodeid().eq.0) then 321 write(6,1) evib, stride, nircopt, inhess, mswg 322 1 format( 323 $ ' energy decrease (evib) = ', 1p,d9.1,0p,/, 324 $ ' initial stride (stride) = ', 1p,d9.1,0p,/, 325 $ ' maximum number of steps (nircopt) = ', i4,/, 326 $ ' initial hessian option (inhess) = ', i4,/, 327 $ ' mass weight coordinates (mswg) = ', l4) 328 write(6,9994) 329 9994 format(/,10x,36('-'), 330 1 /,10x,'Gonzalez & Schlegel IRC Optimization', 331 2 /,10x,36('-'),/) 332c 333 call util_flush(6) 334 end if 335 endif 336c 337c Nvar is the no. of variables in the optimization 338c 339c If we are optimizing the unit cell parameters then we pretend 340c there there are 3 more atoms which will parameterize the 341c unit cell. 342c 343 nat_real = nat 344 ncart = 3*nat 345 nvar = ncart 346 call gsopt_cart_pmat(rtdb, geom) 347c 348 energy = 0d0 349 energyref = 0d0 350 alpha = 1d0 351 gmax = 0d0 352 grms = 0d0 353 smax = 0d0 354 srms = 0d0 355 xmax = 0d0 356 xrms = 0d0 357 call dfill(max_nvar, 0d0, ds, 1) 358 call dfill(max_nvar, 0d0,dsp, 1) 359 call dfill(max_nvar, 0d0, gx, 1) 360 call dfill(max_nvar, 0d0, gq, 1) 361 call dfill(max_nvar, 0d0, g, 1) 362 call dfill(max_nvar, 0d0, oldgra, 1) 363 call dfill(max_nvar, 0d0, sp, 1) 364c 365 if (.not. gsopt_geom_cart_coords_get(geom, sp)) 366 $ call errquit('mepgs: geom?',0, GEOM_ERR) 367c 368 end 369CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCccc 370CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCccc 371 subroutine mepgs_proj_grad(pgref) 372 implicit none 373#include "errquit.fh" 374#include "nwc_const.fh" 375#include "cgsopt.fh" 376#include "cmepgs.fh" 377#include "mafdecls.fh" 378 double precision pgref(nvar) ! returns projected gradient 379c 380c Nothing else is changed. 381c 382 integer l_pmat, k_pmat, l_work, k_work 383c 384 if (.not. ma_push_get(mt_dbl, nvar**2, 'work', 385 $ l_work, k_work)) call errquit 386 $ ('mepgs_proj_h_g: memory for pmat',nvar**2, MA_ERR) 387 if (.not. ma_push_get(mt_dbl, nvar**2, 'pmat', 388 $ l_pmat, k_pmat)) call errquit 389 $ ('mepgs_proj_h_g: memory for work',nvar**2, MA_ERR) 390c 391 call geom_hnd_get_data('p',dbl_mb(k_pmat), nvar**2) 392 if (.not. ma_verify_allocator_stuff()) 393 $ call errquit('freddy',0, MA_ERR) 394c 395c PG 396c 397 call dgemv('n',nvar, nvar, 1d0, dbl_mb(k_pmat), nvar, 398 $ g, 1, 0d0, pgref, 1) 399c 400 if (.not. ma_chop_stack(l_work)) call errquit 401 $ ('mepgs_p_h_g:ma?',0, MA_ERR) 402c 403 end 404CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc 405CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc 406 subroutine mepgs_cent(rtdb, geom, geoma, pgref, sfactor, string) 407 implicit none 408#include "errquit.fh" 409#include "nwc_const.fh" 410#include "cgsopt.fh" 411#include "cmepgs.fh" 412#include "geom.fh" 413#include "rtdb.fh" 414#include "util.fh" 415#include "mafdecls.fh" 416#include "global.fh" 417 integer rtdb, geom, geoma 418 double precision pgref(nvar) 419 double precision sfactor 420 character*(*) string 421c 422c Update the geometry in cent and in the database 423c 'center' by taking the step 424c Update the geometry in geom and in the database 425c 'geometry' by taking the step 426c 427 double precision pgnorm 428 double precision xold(max_cart), xnew(max_cart) 429 integer i, iat, l_bi, k_bi 430 logical gsopt_geom_cart_coords_get 431 logical gsopt_geom_cart_coords_set 432 logical ophigh 433c 434 ophigh = util_print('high', print_high) 435CCCCCCCCCCCCCCCCCCCCCC 436CCCC MASS CCCC 437CCCCCCCCCCCCCCCCCCCCCC 438 if (mswg) call mwgrad(pgref, nvar, .true.) 439CCCCCCCCCCCCCCCCCCCCCC 440CCCC MASS CCCC 441CCCCCCCCCCCCCCCCCCCCCC 442c 443c 444 pgnorm = sqrt(ddot(nvar, pgref, 1, pgref, 1)) 445 call dcopy(nvar, pgref, 1, ds, 1) 446 call dscal(nvar, sfactor*stride/pgnorm, ds, 1) 447c 448c 449CCCCCCCCCCCCCCCCCCCCCC 450CCCC MASS CCCC 451CCCCCCCCCCCCCCCCCCCCCC 452 if (mswg) then 453 call mwgrad(pgref, nvar, .false.) 454 call mwcoord(ds, nvar, .false.) 455 end if 456CCCCCCCCCCCCCCCCCCCCCC 457CCCC MASS CCCC 458CCCCCCCCCCCCCCCCCCCCCC 459c 460c enforce symmetry 461c 462 call gsopt_symmetrize_step(geom) 463c 464c enforce frozen atoms in cartesians 465c 466 if (ga_nodeid().eq.0.and.ophigh) 467 $ write(6,*) 'Zeroing constrained gradient' 468 if ((.not. zcoord) .and. (nactive .ne. nat_real)) then 469 do iat = 1, nat 470 if (.not. oactive(iat)) then 471 do i = 1, 3 472 ds((iat-1)*3+i) = 0.0 473 end do 474 end if 475 end do 476 end if 477c 478 call ga_brdcst(1,ds,8*nvar,0) 479c 480c Get original coordinates 481c 482 if (.not. gsopt_geom_cart_coords_get(geom, xold)) 483 $ call errquit('mepgs_energy_step: coordinates?',geom, 484 & GEOM_ERR) 485c 486 call dcopy(ncart, ds, 1, xnew, 1) 487 call sym_grad_symmetrize(geom, xnew) 488c 489 call daxpy(ncart, 1.0d00, xold, 1, xnew, 1) 490c FRACTIONAL? 491 if (.not. gsopt_geom_cart_coords_set(geoma, xnew)) 492 $ call errquit('mepgs_energy_step: coordinates?', string, 493 & GEOM_ERR) 494c 495 if (.not. geom_rtdb_store(rtdb, geoma, string)) 496 $ call errquit('mepgs_energy_step: grs?',geom, RTDB_ERR) 497c 498 end 499CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 500CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 501 double precision function mepgs_cosang(avec,bvec,angle) 502 implicit none 503#include "errquit.fh" 504#include "nwc_const.fh" 505#include "cgsopt.fh" 506#include "cmepgs.fh" 507#include "util.fh" 508#include "mafdecls.fh" 509c 510 logical angle 511 double precision avec(nvar), bvec(nvar) 512c 513 double precision ctheta, factor(2) 514c 515 mepgs_cosang = 0.0 516c 517 factor(1) = ddot(nvar, avec, 1, bvec, 1) 518c 519 factor(2) = ddot(nvar, avec, 1, avec, 1)* 520 $ ddot(nvar, bvec, 1, bvec, 1) 521c 522 factor(2) = sqrt(factor(2)) 523c 524 ctheta = factor(1)/factor(2) 525c 526 if (abs(ctheta).gt.1.0) ctheta = dsign(1.d0,ctheta) 527c 528 if (angle) then 529 mepgs_cosang = acos(ctheta) 530 else 531 mepgs_cosang = ctheta 532 end if 533c 534 end 535CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 536CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 537 subroutine mepgs_hessian_update() 538 implicit none 539#include "errquit.fh" 540#include "nwc_const.fh" 541#include "cmepgs.fh" 542#include "cgsopt.fh" 543#include "util.fh" 544#include "mafdecls.fh" 545c 546c Update the current Hessian in the optimization variables using 547c . oldgra() - the gradient at the previous point 548c . g() - the gradient at the current point 549c . ds() - the previous search direction 550c . alpha - the step in the previous search direction 551c 552c Only the Hessian is modified. 553c 554 double precision hds(max_nvar) 555 double precision dsds, dshds, dsdg 556 integer l_hess, k_hess, i, j 557 integer ind 558 ind(i,j) = k_hess + i + (j-1)*nvar - 1 559c 560 if (.not. ma_push_get(mt_dbl, nvar**2, 'hess', 561 $ l_hess, k_hess)) call errquit 562 $ ('mepgs_hessian_update: memory for hessian',nvar**2, 563 & GEOM_ERR) 564 call geom_hnd_get_data('irc.hess',dbl_mb(k_hess), nvar**2) 565c 566c Form bits and pieces that are needed 567c 568 call dgemv('n',nvar,nvar,1d0,dbl_mb(k_hess),nvar, 569 $ ds,1,0d0,hds,1) 570c 571 dshds = ddot(nvar, ds, 1, hds, 1) 572 dsds = ddot(nvar, ds, 1, ds, 1) 573 dsdg = 0d0 574 do i = 1, nvar 575 dsdg = dsdg + ds(i)*(g(i) - oldgra(i)) 576 enddo 577c 578c ----- -bfgs- update ----- 579c 580 if(dsdg.gt.1d-14) then 581 do i=1,nvar 582 do j=1,nvar 583 dbl_mb(ind(i,j))=dbl_mb(ind(i,j)) 584 $ + (g(i)-oldgra(i))*(g(j)-oldgra(j))/dsdg 585 1 - hds(i)* hds(j)/dshds 586 enddo 587 enddo 588 endif 589c 590 call geom_hnd_put_data('irc.hess',dbl_mb(k_hess), nvar**2) 591CJMC to be able to start with this hessian OPT 592 call geom_hnd_put_data('gsopt.hess',dbl_mb(k_hess), nvar**2) 593CJMC to be able to start with this hessian OPT 594 if (.not. ma_pop_stack(l_hess)) call errquit 595 $ ('mepgs_hessian_update: ma?',0, MA_ERR) 596c 597 end 598CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 599CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 600 subroutine mepgs_hss_init(rtdb,geom) 601 implicit none 602#include "errquit.fh" 603#include "mafdecls.fh" 604#include "global.fh" 605#include "geom.fh" 606#include "rtdb.fh" 607#include "nwc_const.fh" 608#include "cgsopt.fh" 609#include "cmepgs.fh" 610c 611 integer rtdb, geom 612 613c 614 double precision zero 615 parameter (zero=0.0d+00) 616 integer mxatom, mxcart, mxzmat, mxcoor 617 parameter (mxatom=nw_max_atom) 618 parameter (mxcart=3*mxatom) 619 parameter (mxzmat=nw_max_zmat) 620 parameter (mxcoor=nw_max_coor) 621c 622c These commons are used in the internal coordinate guess 623c 624 integer nuc 625 COMMON/HND_MOLNUC/NUC(MXATOM) 626 double precision c, zan 627 integer natom 628 common/hnd_molxyz/c(3,mxatom),zan(mxatom),natom 629 integer nnzmat, nnzvar, nnvar 630 common/hnd_zmtpar/nnzmat,nnzvar,nnvar 631 double precision hscale, ascale, bscale, tscale, amat(3,3) 632c 633 integer l_hess, k_hess, l_zmat, k_zmat, l_izmat, k_izmat, i, j 634 integer l_c, k_c, l_t, k_t, iat 635 logical old_hessian 636 character*16 atom_tags(mxatom) 637c 638 nnzmat = nzmat 639 nnzvar = nzvar 640 nnvar = nzvar 641 if (.not. geom_ncent(geom,natom)) 642 1 call errquit('hnd_opt: geom_ncent?',911, GEOM_ERR) 643c 644 if (.not. ma_push_get(mt_dbl, nvar**2, 'hessian', 645 $ l_hess, k_hess)) call errquit 646 $ ('mepgs_init_hess: failed allocating hessian',nvar**2, 647 & MA_ERR) 648c 649 old_hessian=.false. 650 call gsopt_check_hess(nvar, old_hessian) 651 if (oprint) write(6,*) 652 if (old_hessian) then 653 call mepgs_hess_cart_guess() 654 if (oprint) write(6,*) 655 $ ' Using Cartesian Hessian from previous frequency', 656 $ ' calculation' 657 else 658 if (oprint) write(6,*) 'Not restart Hessian? ' 659 endif 660c 661c Apply constants, constraints and overall scaling 662c 663 if (.not. ma_push_get(mt_dbl, nvar**2, 'hessian', 664 $ l_c, k_c)) call errquit 665 $ ('mepgs_init_hess: failed allocating hessian',nvar**2, 666 & MA_ERR) 667 if (.not. ma_push_get(mt_dbl, nvar**2, 'hessian', 668 $ l_t, k_t)) call errquit 669 $ ('mepgs_init_hess: failed allocating hessian',nvar**2, 670 & MA_ERR) 671c 672 call geom_hnd_get_data('irc.hess', dbl_mb(k_hess), nvar*nvar) 673c 674c Used to use c here ... now use p 675c 676 call geom_hnd_get_data('p',dbl_mb(k_c), nvar*nvar) 677c 678 if (odebug) then 679 write(6,*) ' Initial Hessian before P' 680 call output(dbl_mb(k_hess),1,nvar,1,nvar,nvar,nvar,1) 681 endif 682 call dgemm('n','n',nvar,nvar,nvar,1d0,dbl_mb(k_c),nvar, 683 $ dbl_mb(k_hess),nvar,0d0,dbl_mb(k_t),nvar) 684 call dgemm('n','t',nvar,nvar,nvar,1d0,dbl_mb(k_t),nvar, 685 $ dbl_mb(k_c),nvar,0d0,dbl_mb(k_hess),nvar) 686 if (odebug) then 687 write(6,*) ' Initial Hessian after P' 688 call output(dbl_mb(k_hess),1,nvar,1,nvar,nvar,nvar,1) 689 endif 690c 691 if (nactive .ne. nat_real) then 692c 693c We are in cartesian coordinates and some have been frozen. 694c Since there is no redundancy or coupling we just need 695c to make sure that the initial Hessian does not couple 696c frozen with unfrozen variables and we are OK. 697c 698 do iat = 1, nat 699 if (.not. oactive(iat)) then 700 do i = 1+(iat-1)*3, iat*3 701 do j = 1, nvar 702 dbl_mb(k_hess+j-1+(i-1)*nvar) = 0d0 703 dbl_mb(k_hess+i-1+(j-1)*nvar) = 0d0 704 enddo 705 dbl_mb(k_hess+i-1+(i-1)*nvar) = 1d0 706 enddo 707 endif 708 enddo 709 endif 710c 711 if (.not. rtdb_get(rtdb,'gsopt:hscale',mt_dbl,1,hscale)) 712 $ hscale = 1d0 713 call dscal(nvar*nvar, hscale, dbl_mb(k_hess), 1) 714 if (oprint .and. hscale.ne.1d0) then 715 if (ga_nodeid().eq.0) write(6,78) hscale 716 end if 717 78 format(' Scaling initial hessian by ',f6.2) 718c 719 call geom_hnd_put_data('irc.hess',dbl_mb(k_hess), nvar*nvar) 720CJMC prepare OPT start hessian 721 call geom_hnd_put_data('gsopt.hess',dbl_mb(k_hess), nvar*nvar) 722CJMC prepare OPT start hessian 723c 724 if (.not. ma_chop_stack(l_hess)) call errquit 725 $ ('mepgs_init_hess ma corrupt',0, MA_ERR) 726c 727 END 728CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 729CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 730 subroutine mepgs_hess_cart_guess() 731 implicit none 732#include "errquit.fh" 733#include "mafdecls.fh" 734#include "global.fh" 735#include "nwc_const.fh" 736#include "cmepgs.fh" 737#include "cgsopt.fh" 738#include "inp.fh" 739c 740c Read in cartesian Hessian and transform it as necessary 741c to internal coordinates (neglecting the component due to 742c the derivative) and writing the result to the hessian file. 743c 744c Reads file in vib_vib format using vib_vib filename default 745c Note the default filename is set in task_freq 746c filenames must be made identical. 747c 748c Format of vib file is ascii lower triangular elements only. 749c 750 integer h_unit 751 parameter (h_unit=47) 752 character*255 fname 753 double precision x 754 integer i,j 755 integer l_bi, k_bi, l_hc, k_hc, l_hq, k_hq 756c 757 if (.not. ma_push_get(mt_dbl, ncart*nvar, 'binv', 758 $ l_bi, k_bi)) call errquit 759 $ ('mepgs_hess_cart_guess: ma?', ncart*nvar, MA_ERR) 760c 761 if (.not. ma_push_get(mt_dbl, max(ncart**2,nvar**2), 'hcart', 762 $ l_hc, k_hc)) call errquit 763 $ ('mepgs_hess_cart_guess: ma?', ncart**2, MA_ERR) 764c 765 if (.not. ma_push_get(mt_dbl, max(ncart**2,nvar**2), 'hcart2', 766 $ l_hq, k_hq)) call errquit 767 $ ('mepgs_hess_cart_guess: ma?', nvar**2, MA_ERR) 768c 769 if (ga_nodeid().eq.0) then 770 call util_file_name('hess',.false.,.false.,fname) 771 open(unit=h_unit,file=fname,form='formatted',status='unknown', 772 $ err=99990,access='sequential') 773 rewind h_unit 774 do i = 1,ncart 775 do j = 1,i 776 read(h_unit,10000,err=99992,end=99992) x 777 dbl_mb(k_hc+(i-1)*ncart+(j-1)) = x 778 dbl_mb(k_hc+(j-1)*ncart+(i-1)) = x 779 enddo 780 enddo 781 close(unit=h_unit,status='keep') 782 endif 783 call ga_brdcst(1,dbl_mb(k_hc),8*ncart**2,0) 784c 785 call geom_hnd_get_data('b^-1', dbl_mb(k_bi), nvar*ncart) 786 call dgemm('n', 'n', ncart, nvar, ncart, 1d0, dbl_mb(k_hc), ncart, 787 $ dbl_mb(k_bi), ncart, 0d0, dbl_mb(k_hq), ncart) 788 call dgemm('t', 'n', nvar, nvar, ncart, 1d0, dbl_mb(k_bi), ncart, 789 $ dbl_mb(k_hq), ncart, 0d0, dbl_mb(k_hc), nvar) 790c 791 do i = 1,nvar 792 do j = 1,i 793 x = (dbl_mb(k_hc+(i-1)*nvar+(j-1)) + 794 $ dbl_mb(k_hc+(j-1)*nvar+(i-1))) * 0.5d0 795 dbl_mb(k_hc+(i-1)*nvar+(j-1)) = x 796 dbl_mb(k_hc+(j-1)*nvar+(i-1)) = x 797 enddo 798 enddo 799c 800CCCCCCCCCCCCCCCCCCCCCC 801CCCC MASS CCCC 802CCCCCCCCCCCCCCCCCCCCCC 803 if (mswg) call geom_hnd_get_data('irc.hess',dbl_mb(k_hc), nvar**2) 804CCCCCCCCCCCCCCCCCCCCCC 805CCCC MASS CCCC 806CCCCCCCCCCCCCCCCCCCCCC 807 call geom_hnd_put_data('irc.hess',dbl_mb(k_hc), nvar**2) 808c 809 if (.not. ma_chop_stack(l_bi)) 810 $ call errquit('mepgs_hess_cart_guess: ma corrupt?',0, MA_ERR) 811c 812 return 81310000 format(f30.15) 81499990 write(6,*)' could not open <',fname(1:inp_strlen(fname)), 815 $ '> as unknown file' 816 call errquit('mepgs_hess_cart: fatal error', 911, GEOM_ERR) 81799991 write(6,*)' could not open <',fname(1:inp_strlen(fname)), 818 $ '> as new file' 819 call errquit('mepgs_hess_cart: fatal error', 911, GEOM_ERR) 82099992 write(6,*)' error in reading <',fname(1:inp_strlen(fname)), 821 $ '> as hessian file' 822 call errquit('mepgs_hess_cart: fatal error', 911, GEOM_ERR) 823 end 824CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 825CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 826 subroutine mepgs_path(geom, openfile, closefile) 827 implicit none 828#include "errquit.fh" 829#include "nwc_const.fh" 830#include "cgsopt.fh" 831#include "cmepgs.fh" 832#include "global.fh" 833#include "util.fh" 834#include "inp.fh" 835 integer geom 836 logical openfile, closefile 837 character*255 filename, dir 838 logical mol_geom_print_xyz 839 external mol_geom_print_xyz 840c 841c Print a trajectory file 842c 843 if (ga_nodeid().eq.0 .and. xyz.ne.' ') then 844 dir = ' ' 845 filename = ' ' 846 call util_directory_name(dir, .false., 0) 847CJMC 848 if (forward) then 849 write(filename,12) dir(1:inp_strlen(dir)), 850 $ xyz(1:inp_strlen(xyz)) 851 12 format(a,'/',a,'.fxyz') 852 else if (backward) then 853 write(filename,13) dir(1:inp_strlen(dir)), 854 $ xyz(1:inp_strlen(xyz)) 855 13 format(a,'/',a,'.bxyz') 856 end if 857CJMC 858 if (openfile) then 859 open(88,file=filename,form='formatted',status='unknown', 860 $ access='sequential',err=133) 861 rewind(88) 862 end if 863c 864 if (.not. mol_geom_print_xyz(geom, 88, energy)) 865 $ call errquit('mepgs_path: mol_geom_print_xyz?',0, GEOM_ERR) 866 call util_flush(88) 867c 868 if (closefile) close(88,status='keep',err=133) 869 end if 870c 871 return 872 133 call errquit('mepgs_path: error open/close xyz file',0, GEOM_ERR) 873c 874 end 875CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc 876 subroutine mepgs_gra(rtdb, geom) 877 implicit none 878#include "errquit.fh" 879#include "global.fh" 880#include "mafdecls.fh" 881#include "stdio.fh" 882#include "util.fh" 883#include "nwc_const.fh" 884#include "cgsopt.fh" 885#include "cmepgs.fh" 886#include "rtdb.fh" 887#include "geom.fh" 888c 889 integer rtdb, geom 890c 891 integer iat, ixyz 892 logical ophigh 893 logical task_gradient 894 external task_gradient 895c 896 ophigh = util_print('high', print_high) 897c 898 if (.not. task_gradient(rtdb)) 899 $ call errquit('mepgs: task_gradient failed',0, GEOM_ERR) 900 call gsopt_get_grad(rtdb, geom) ! Into gx 901c 902c Zero the gradient associated with atoms frozen in cartesians 903c 904 if (ga_nodeid().eq.0.and.ophigh) 905 $ write(6,*) 'Zeroing constrained gradient' 906 if ((.not. zcoord) .and. (nactive .ne. nat_real)) then 907 do iat = 1, nat_real 908 if (.not. oactive(iat)) then 909 do ixyz = 1, 3 910 gx((iat-1)*3+ixyz) = 0.0 911 end do 912 end if 913 end do 914 end if 915c 916 end 917CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 918 subroutine mepgs_chk(geom, irccyc) 919 implicit none 920#include "errquit.fh" 921#include "global.fh" 922#include "mafdecls.fh" 923#include "stdio.fh" 924#include "util.fh" 925#include "nwc_const.fh" 926#include "cgsopt.fh" 927#include "cmepgs.fh" 928#include "geom.fh" 929c 930 integer geom, irccyc 931c 932 integer ivar 933 double precision abvec(max_nvar) 934 double precision bcvec(max_nvar) 935 double precision newx(max_nvar) 936 double precision ogrms 937 double precision graang, cosang, stpang, rmsdif 938 double precision eps 939 parameter (eps = 1d-4) 940 double precision echange 941c 942 character*1 mark 943c 944 double precision mepgs_cosang 945 logical gsopt_geom_cart_coords_get 946c 947 948 if (.not. gsopt_geom_cart_coords_get(geom, newx)) 949 $ call errquit('mepgs: geom?',0, geom_err) 950c 951c 952CCCCCCCCCCCCCCCCCCCCCC 953CCCC MASS CCCC 954CCCCCCCCCCCCCCCCCCCCCC 955 if (mswg) then 956 call mwcoord( ds, nvar, .true.) 957 call mwcoord( newx, nvar, .true.) 958 call mwcoord(center, nvar, .true.) 959 call mwcoord(oldgeo, nvar, .true.) 960 call mwgrad(oldgra, nvar, .true.) 961 call mwgrad( g, nvar, .true.) 962 end if 963CCCCCCCCCCCCCCCCCCCCCC 964CCCC MASS CCCC 965CCCCCCCCCCCCCCCCCCCCCC 966c 967c 968 call dcopy(nvar, center, 1, abvec, 1) 969 call daxpy(nvar, -1.0d0, oldgeo, 1, abvec, 1) 970c 971 call dcopy(nvar, newx, 1, bcvec, 1) 972 call daxpy(nvar, -1.0d0, center, 1, bcvec, 1) 973c 974 graang = mepgs_cosang(g, oldgra, .true.) 975 cosang = mepgs_cosang(g, oldgra, .false.) 976 stpang = mepgs_cosang(abvec, bcvec, .true.) 977c 978 grms = 0d0 979 srms = 0d0 980 gmax = 0d0 981 smax = 0d0 982c 983 do ivar = 1, nvar 984 ogrms = ogrms + oldgra(ivar)*oldgra(ivar) 985 grms = grms + g(ivar)*g(ivar) 986 gmax = max(gmax, abs(g(ivar))) 987 enddo 988 ogrms = sqrt(ogrms/dble(nvar)) 989 grms = sqrt(grms/dble(nvar)) 990c 991 rmsdif = grms - ogrms 992 if (rmsdif.lt.0d0) flip = .true. 993 if (abs(rmsdif).lt.eps) rmsdif = 0d0 994 if (irccyc.eq.1) flip = .false. 995c 996 svalue = graang*0.5d0*stride*(1d0/tan(graang/2.0d0)) 997c 998 echange = energy - energyref 999c 1000 if (flip .and. (grms.le.grms_tol) .and. (gmax.le.gmax_tol) 1001 $ .and. (abs(echange).lt.eprec)) ircdone = .true. 1002c 1003c Assume step was good 1004c 1005 redogs = .false. 1006c 1007 if (flip) then 1008 if((rmsdif.gt.0d0).or.(echange.gt.0d0)) redogs=.true. 1009 else 1010 if((echange.gt.0d0).and.(abs(echange).gt.eps)) redogs=.true. 1011 end if 1012c 1013 if ((irccyc.ne.1).and.(cosang.le.-0.7d0)) redogs = .true. 1014c 1015 if (.not. redogs) then 1016 stotal = stotal + svalue 1017 mark = '&' 1018 if (ga_nodeid().eq.0) write(6,1) mark, mark 1019 mark = '&' 1020 if (ga_nodeid().eq.0) 1021 & write(6,2) mark, irccyc-1, energy, svalue, stotal 1022 1 format( 1023 $ /,a1,' Point Energy svalue stotal ', 1024 $ /,a1,' ------- ------------ -------- --------') 1025 2 format( 1026 $ a1,i5,f17.8,2f9.5,/) 1027 end if 1028c 1029c 1030CCCCCCCCCCCCCCCCCCCCCC 1031CCCC MASS CCCC 1032CCCCCCCCCCCCCCCCCCCCCC 1033 if (mswg) then 1034 call mwcoord( ds, nvar, .false.) 1035 call mwcoord( newx, nvar, .false.) 1036 call mwcoord(center, nvar, .false.) 1037 call mwcoord(oldgeo, nvar, .false.) 1038 call mwgrad(oldgra, nvar, .false.) 1039 call mwgrad( g, nvar, .false.) 1040 end if 1041CCCCCCCCCCCCCCCCCCCCCC 1042CCCC MASS CCCC 1043CCCCCCCCCCCCCCCCCCCCCC 1044c 1045c 1046 end 1047CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc 1048 subroutine mepgs_dist(geom) 1049 implicit none 1050#include "errquit.fh" 1051#include "global.fh" 1052#include "mafdecls.fh" 1053#include "stdio.fh" 1054#include "util.fh" 1055#include "nwc_const.fh" 1056#include "cgsopt.fh" 1057#include "cmepgs.fh" 1058#include "geom.fh" 1059c 1060 integer geom 1061c 1062 integer ivar 1063 double precision newx(max_nvar) 1064 logical gsopt_geom_cart_coords_get 1065c 1066 if (.not. gsopt_geom_cart_coords_get(geom, newx)) 1067 $ call errquit('mepgs: geom?',0, geom_err) 1068c 1069c 1070CCCCCCCCCCCCCCCCCCCCCC 1071CCCC MASS CCCC 1072CCCCCCCCCCCCCCCCCCCCCC 1073 if (mswg) then 1074 call mwcoord(center, nvar, .true.) 1075 call mwcoord( newx, nvar, .true.) 1076 end if 1077CCCCCCCCCCCCCCCCCCCCCC 1078CCCC MASS CCCC 1079CCCCCCCCCCCCCCCCCCCCCC 1080c 1081c 1082 distcx = 0d0 1083 do ivar=1, nvar 1084 distcx = distcx + (center(ivar) - newx(ivar))**2 1085 end do 1086c 1087c 1088CCCCCCCCCCCCCCCCCCCCCC 1089CCCC MASS CCCC 1090CCCCCCCCCCCCCCCCCCCCCC 1091 if (mswg) then 1092 call mwcoord(center, nvar, .false.) 1093 call mwcoord( newx, nvar, .false.) 1094 end if 1095CCCCCCCCCCCCCCCCCCCCCC 1096CCCC MASS CCCC 1097CCCCCCCCCCCCCCCCCCCCCC 1098c 1099c 1100 end 1101CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1102 logical function mepgs_opt(rtdb) 1103 implicit none 1104#include "errquit.fh" 1105#include "global.fh" 1106#include "mafdecls.fh" 1107#include "stdio.fh" 1108#include "util.fh" 1109#include "nwc_const.fh" 1110#include "cgsopt.fh" 1111#include "cmepgs.fh" 1112#include "rtdb.fh" 1113#include "geom.fh" 1114c 1115 integer rtdb 1116c 1117 integer geom 1118 integer reference 1119 integer imepgs 1120 logical ignore 1121 double precision mepgs_cosang 1122 logical status 1123 double precision start ! Tracks time used in last step 1124cjmc 1125 double precision minst ! Tracks time used in last step 1126 parameter(minst = 0.002d0) 1127 double precision pgref(max_nvar) 1128 logical gsopt_geom_cart_coords_get 1129 logical gsopt_geom_cart_coords_set 1130 logical gsopt 1131cjmc 1132 integer required ! Time required 1133 logical ophigh 1134 logical task_gradient, task_energy 1135 external task_gradient, task_energy 1136c 1137c **** Initialization **** 1138c 1139 if (.not. geom_create(geom, 'geometry')) 1140 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 1141 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 1142 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 1143 if (.not. geom_ncent(geom,nat)) 1144 $ call errquit('hnd_opt: natoms?',nat, GEOM_ERR) 1145 call grad_active_atoms(rtdb, nat, oactive, nactive) 1146 if (.not. geom_systype_get(geom, isystype)) 1147 $ call errquit('mepgs: systype?',0, GEOM_ERR) 1148 call gsopt_get_grad(rtdb, geom) ! Into gx 1149c 1150c **** Gonzalez & Schlegel Iterative loop **** 1151c 1152 do imepgs = 1, nircopt+1 1153c 1154c **** Iteration exceded **** 1155c 1156 if (imepgs.gt.nircopt) goto 200 1157c 1158 start = util_wallsec() 1159 if (oprint.and.ga_nodeid().eq.0) write(6,1) imepgs-1 1160 1 format(/,10x,11('-'),/,10x,'GS Step',i4,/,10x,11('-')) 1161 if ((ga_nodeid() .eq. 0) .and. 1162 $ util_print('geometry',print_default)) then 1163 if (.not. geom_print(geom)) call errquit('mepgs: geom?',0, 1164 $ GEOM_ERR) 1165 endif 1166c 1167c Save old energy, gradient and coordinates 1168c 1169 energyref = energy 1170 call dcopy(ncart, gx, 1, g, 1) ! g() set to gx() 1171 call dcopy(nvar, g, 1, oldgra, 1) 1172 call dcopy(nvar, g, 1, gp, 1) 1173 if (.not. gsopt_geom_cart_coords_get(geom, oldgeo)) 1174 $ call errquit('mepgs: geom?',0, geom_err) 1175c 1176c Project gradient 1177c 1178 call mepgs_proj_grad(pgref) 1179c 1180 1000 continue 1181c 1182c **** Allocate and prepare field for reference (center) **** 1183c 1184 if (.not. geom_create(reference, 'center')) 1185 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 1186 if (.not. geom_rtdb_load(rtdb, reference, 'geometry')) 1187 & call errquit('hnd_opt: no initial geometry ',911,RTDB_ERR) 1188c 1189c **** Calculate center, half stride **** 1190c 1191 call mepgs_cent(rtdb, geom, reference, pgref, -0.5d0,'center') 1192c 1193c **** Obtain center coordinates **** 1194c 1195 if (.not. gsopt_geom_cart_coords_get(reference, center)) 1196 $ call errquit('mepgs: geom?',0, geom_err) 1197c 1198c **** Deallocate reference (center) **** 1199c 1200 if (.not. geom_destroy(reference)) 1201 $ call errquit('mepgs:reference corrupt',0, GEOM_ERR) 1202c 1203c **** Calculate x^(l), complete stride **** 1204c 1205 call mepgs_cent(rtdb, geom, geom, pgref, -1.0d0, 'geometry') 1206c 1207c **** Energy and Gradient **** 1208c 1209 call mepgs_gra(rtdb, geom) 1210 call dcopy(ncart, gx, 1, g, 1) ! g() set to gx() 1211c 1212c **** Reject proposed point if necessary **** 1213c **** Based on gradient angle **** 1214c 1215 if (flip) then 1216c 1217c 1218CCCCCCCCCCCCCCCCCCCCCC 1219CCCC MASS CCCC 1220CCCCCCCCCCCCCCCCCCCCCC 1221 if (mswg) then 1222 call mwgrad( g, nvar, .true.) 1223 call mwgrad(oldgra, nvar, .true.) 1224 end if 1225CCCCCCCCCCCCCCCCCCCCCC 1226CCCC MASS CCCC 1227CCCCCCCCCCCCCCCCCCCCCC 1228c 1229c 1230 if (mepgs_cosang(g, oldgra, .false.).le.-0.7d0) then 1231 if (oprint) write(6,*) "Rejecting the proposed Point" 1232 if (stride.eq.minst) then 1233 ircdone = .true. 1234 goto 100 1235 end if 1236 stride = stride/2d0 1237 stride = max(stride,minst) 1238 call dcopy(nvar, oldgra, 1, g, 1) 1239c 1240c 1241CCCCCCCCCCCCCCCCCCCCCC 1242CCCC MASS CCCC 1243CCCCCCCCCCCCCCCCCCCCCC 1244 if (mswg) then 1245 call mwgrad( g, nvar, .false.) 1246 call mwgrad(oldgra, nvar, .false.) 1247 end if 1248CCCCCCCCCCCCCCCCCCCCCC 1249CCCC MASS CCCC 1250CCCCCCCCCCCCCCCCCCCCCC 1251c 1252c 1253 if (.not. gsopt_geom_cart_coords_set(geom, oldgeo)) 1254 $ call errquit('gsopt: geom?',0, geom_err) ! reload previous cart 1255 if (.not. geom_rtdb_store(rtdb, geom, 'geometry')) 1256 $ call errquit('gsopt: grs?',geom, RTDB_ERR) 1257 goto 1000 1258 end if 1259c 1260c 1261CCCCCCCCCCCCCCCCCCCCCC 1262CCCC MASS CCCC 1263CCCCCCCCCCCCCCCCCCCCCC 1264 if (mswg) then 1265 call mwgrad( g, nvar, .false.) 1266 call mwgrad(oldgra, nvar, .false.) 1267 end if 1268CCCCCCCCCCCCCCCCCCCCCC 1269CCCC MASS CCCC 1270CCCCCCCCCCCCCCCCCCCCCC 1271c 1272c 1273 end if 1274c 1275c **** Update the OPT hessian **** 1276c 1277 call dcopy(nvar, oldgeo, 1, sp, 1) 1278 call gsopt_compute_actual_step(geom) 1279c 1280c 1281CCCCCCCCCCCCCCCCCCCCCC 1282CCCC MASS CCCC 1283CCCCCCCCCCCCCCCCCCCCCC 1284 if (mswg) then 1285 call mwcoord( ds, nvar, .true.) 1286 call mwgrad( g, nvar, .true.) 1287 call mwgrad( gp, nvar, .true.) 1288 end if 1289CCCCCCCCCCCCCCCCCCCCCC 1290CCCC MASS CCCC 1291CCCCCCCCCCCCCCCCCCCCCC 1292c 1293c 1294 call gsopt_hessian_update() 1295c 1296c 1297CCCCCCCCCCCCCCCCCCCCCC 1298CCCC MASS CCCC 1299CCCCCCCCCCCCCCCCCCCCCC 1300 if (mswg) then 1301 call mwcoord( ds, nvar, .false.) 1302 call mwgrad( g, nvar, .false.) 1303 call mwgrad( gp, nvar, .false.) 1304 end if 1305CCCCCCCCCCCCCCCCCCCCCC 1306CCCC MASS CCCC 1307CCCCCCCCCCCCCCCCCCCCCC 1308c 1309c 1310c **** Compute initial distance **** 1311c 1312 call mepgs_dist(geom) 1313 ctrust2 = distcx 1314c 1315c **** Deallocate geom **** 1316c 1317 if (.not.geom_destroy(geom)) 1318 & call errquit('hnd_opt: geom_destroy?', 911, GEOM_ERR) 1319c 1320c **** Perform constrained optimization **** 1321c 1322 if (.not. gsopt(rtdb)) 1323 $ call errquit('mepgs: gsopt failed',0, GEOM_ERR) 1324c 1325c *** Reallocate geom info **** 1326c 1327 if (.not. geom_create(geom, 'geometry')) 1328 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 1329 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 1330 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 1331 if (.not. geom_ncent(geom,nat)) 1332 $ call errquit('hnd_opt: natoms?',nat, GEOM_ERR) 1333 call grad_active_atoms(rtdb, nat, oactive, nactive) 1334 if (.not. geom_systype_get(geom, isystype)) 1335 $ call errquit('mepgs: systype?',0, GEOM_ERR) 1336c 1337c **** Compute final distance **** 1338c 1339 call mepgs_dist(geom) 1340 if (abs(distcx - ctrust2) .gt. 1d-4) 1341 $ call errquit('mepgs: norm not preserved', 911, GEOM_ERR) 1342c 1343c **** obtain current gradient **** 1344c 1345 call gsopt_get_grad(rtdb, geom) ! Into gx 1346 call dcopy(ncart, gx, 1, g, 1) ! g() set to gx() 1347c 1348c **** Compute convergence info **** 1349c 1350 call mepgs_chk(geom, imepgs) 1351c 1352c **** reload previous fiels if necessary **** 1353c 1354 if (redogs) then 1355 if (oprint) write(6,*) "Rejecting the optimized Point" 1356 if (stride.eq.minst) then 1357 ircdone = .true. 1358 goto 100 1359 end if 1360 stride = stride/2d0 1361 stride = max(stride,minst) 1362 if (.not. gsopt_geom_cart_coords_set(geom, oldgeo)) 1363 $ call errquit('gsopt: geom?',0, geom_err) ! reload previous cart 1364 if (.not. geom_rtdb_store(rtdb, geom, 'geometry')) 1365 $ call errquit('gsopt: grs?',geom, RTDB_ERR) 1366 goto 1000 1367c 1368 end if 1369c 1370c Check for convergence 1371c 1372 if (ircdone) goto 100 1373c 1374c Update the IRC hessian 1375c 1376 call dcopy(nvar, oldgeo, 1, sp, 1) 1377 call gsopt_compute_actual_step(geom) 1378c 1379c 1380CCCCCCCCCCCCCCCCCCCCCC 1381CCCC MASS CCCC 1382CCCCCCCCCCCCCCCCCCCCCC 1383 if (mswg) then 1384 call mwcoord( ds, nvar, .true.) 1385 call mwgrad( g, nvar, .true.) 1386 call mwgrad(oldgra, nvar, .true.) 1387 end if 1388CCCCCCCCCCCCCCCCCCCCCC 1389CCCC MASS CCCC 1390CCCCCCCCCCCCCCCCCCCCCC 1391c 1392c 1393 call mepgs_hessian_update() 1394c 1395c 1396CCCCCCCCCCCCCCCCCCCCCC 1397CCCC MASS CCCC 1398CCCCCCCCCCCCCCCCCCCCCC 1399 if (mswg) then 1400 call mwcoord( ds, nvar, .false.) 1401 call mwgrad( g, nvar, .false.) 1402 call mwgrad(oldgra, nvar, .false.) 1403 end if 1404CCCCCCCCCCCCCCCCCCCCCC 1405CCCC MASS CCCC 1406CCCCCCCCCCCCCCCCCCCCCC 1407c 1408c 1409c Print a trajectory file 1410c 1411 call mepgs_path(geom, .false., .false.) 1412c 1413c Check time before next iteration 1414c 1415 required = int(1.2d0*(util_wallsec() - start)) + 1 1416 if (.not. util_test_time_remaining(rtdb,required)) goto 200 1417c 1418c 1419 enddo ! End of iterative loop 1420c 1421c **** Failed to converge **** 1422c 1423 200 if (oprint.and.ga_nodeid().eq.0) write(6,201) 1424 201 format(/,1x,63('-')/,1x,'Failed to converge in maximum number', 1425 $ ' of steps or available time'/,1x,63('-')/) 1426 ircdone = .false. 1427c 1428c **** Procedure finished **** 1429c 1430 100 if (ircdone) then 1431 if (oprint.and.ga_nodeid().eq.0) write(6,101) 1432 101 format(/,6x,22('-'),/,6x,'IRC Optimization converged',/, 1433 $ 6x,22('-'),/) 1434 endif 1435c 1436c Print out final info and geometry 1437c 1438 if (ga_nodeid().eq.0 .and. util_print('finish',print_low)) then 1439 call mepgs_path(geom, .false., .true.) 1440 if (.not. geom_print(geom)) call errquit 1441 $ ('hnd_opt_drv: geom_print?',0, GEOM_ERR) 1442c 1443 end if 1444c 1445c **** Save geometry **** 1446c 1447 if (.not. geom_rtdb_store(rtdb, geom, 'geometry')) 1448 $ call errquit('gsopt: grs?',geom, RTDB_ERR) 1449c 1450c Clean up and go home 1451c 1452 if (.not.geom_destroy(geom)) 1453 & call errquit('hnd_opt: geom_destroy?', 911, GEOM_ERR) 1454c 1455 mepgs_opt=ircdone 1456 if (ircdone) then 1457 call ecce_print_module_exit('mepgs', 'ok') 1458 else 1459 call ecce_print_module_exit('mepgs', 'failed') 1460 endif 1461c 1462 call movecs_ecce_print_on() ! Restore MO printing 1463 call util_print_pop 1464c 1465 call ga_sync() 1466c 1467 end 1468CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1469CCCCCCCCC GS optimization CCCCCCCCCC 1470CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1471 logical function gsopt(rtdb) 1472 implicit none 1473#include "errquit.fh" 1474#include "global.fh" 1475#include "mafdecls.fh" 1476#include "stdio.fh" 1477#include "util.fh" 1478#include "nwc_const.fh" 1479#include "cgsopt.fh" 1480#include "cmepgs.fh" 1481#include "rtdb.fh" 1482#include "geom.fh" 1483c 1484 integer rtdb 1485c 1486 integer geom, i, iat 1487 integer istep 1488 logical converged, status 1489 double precision start ! Tracks time used in last step 1490cjmc 1491 double precision oldx(max_cart) 1492 logical gsopt_geom_cart_coords_get 1493 logical gsopt_geom_cart_coords_set 1494cjmc 1495 integer required ! Time required 1496 logical ophigh 1497 logical gsopt_converged, task_gradient 1498 external gsopt_converged, task_gradient 1499c 1500 ophigh = util_print('high', print_high) 1501c 1502c Read input, load /cgsopt/, get geometry 1503c 1504 call gsopt_initialize(rtdb, geom) 1505c 1506c Energy and Gradient 1507c 1508 call gsopt_get_grad(rtdb, geom) ! Into gx 1509 if (.not. rtdb_get(rtdb,'task:energy', mt_dbl, 1, energy)) 1510 $ call errquit('gsopt: could not get energy',0, RTDB_ERR) 1511c 1512c Construct proyector 1513c 1514 call gsopt_cart_pmat(rtdb, geom) 1515 call dcopy(ncart, gx, 1, g, 1) ! g() set to gx() 1516c 1517c Initial hessian 1518c 1519 call gsopt_hss_init(rtdb,geom) 1520c 1521c Iterative loop 1522c 1523 do istep = 1, nptopt+1 ! +1 since first pass thru loop is not a step 1524c 1525c iteration exceded 1526c 1527 if(istep.gt.nptopt) goto 200 1528c 1529 start = util_wallsec() 1530 if (oprint.and. ga_nodeid().eq.0) write(6,1) istep-1 1531 1 format(/,10x,'--------',/,10x,'Step',i4,/,10x,'--------') 1532 if ((ga_nodeid() .eq. 0) .and. 1533 $ util_print('geometry',print_default)) then 1534 if (.not. geom_print(geom)) call errquit('gsopt: geom?',0, 1535 $ GEOM_ERR) 1536 endif 1537c 1538c Compute step/gradient info and print for user 1539c (for current energy & gradient, and the previous alpha*step). 1540c 1541c 1542CCCCCCCCCCCCCCCCCCCCCC 1543CCCC MASS CCCC 1544CCCCCCCCCCCCCCCCCCCCCC 1545 if (mswg) then 1546 call mwcoord(radius, nvar, .true.) 1547 call mwcoord( ds, nvar, .true.) 1548 call mwgrad( g, nvar, .true.) 1549 end if 1550CCCCCCCCCCCCCCCCCCCCCC 1551CCCC MASS CCCC 1552CCCCCCCCCCCCCCCCCCCCCC 1553c 1554c 1555 call gsopt_compute_info() 1556c 1557c 1558CCCCCCCCCCCCCCCCCCCCCC 1559CCCC MASS CCCC 1560CCCCCCCCCCCCCCCCCCCCCC 1561 if (mswg) then 1562 call mwcoord(radius, nvar, .false.) 1563 call mwcoord( ds, nvar, .false.) 1564 call mwgrad( g, nvar, .false.) 1565 end if 1566CCCCCCCCCCCCCCCCCCCCCC 1567CCCC MASS CCCC 1568CCCCCCCCCCCCCCCCCCCCCC 1569c 1570c 1571 call gsopt_print(geom, istep) 1572c 1573c Check for convergence 1574c 1575 if (gsopt_converged(istep)) then 1576 converged = .true. 1577 goto 100 1578 endif 1579c 1580c Save old energy, gradient and coordinates 1581C 1582 energyp = energy ! Used for convergence and step restriction 1583 call dcopy(nvar, g, 1, gp, 1) ! Used for Hessian update 1584 if (.not. gsopt_geom_cart_coords_get(geom, oldx)) 1585 $ call errquit('gsopt: geom?',0, geom_err) 1586c 1587c generate a new search direction 1588c 1589c 1590CCCCCCCCCCCCCCCCCCCCCC 1591CCCC MASS CCCC 1592CCCCCCCCCCCCCCCCCCCCCC 1593 if (mswg) call mwgrad(g, nvar, .true.) 1594CCCCCCCCCCCCCCCCCCCCCC 1595CCCC MASS CCCC 1596CCCCCCCCCCCCCCCCCCCCCC 1597c 1598c 1599 call gsopt_pickstp(rtdb, geom, istep) ! fills in ds() 1600c 1601c 1602CCCCCCCCCCCCCCCCCCCCCC 1603CCCC MASS CCCC 1604CCCCCCCCCCCCCCCCCCCCCC 1605 if (mswg) call mwgrad(g, nvar, .false.) 1606CCCCCCCCCCCCCCCCCCCCCC 1607CCCC MASS CCCC 1608CCCCCCCCCCCCCCCCCCCCCC 1609c 1610c take recommended step. 1611c 1612 call gsopt_take_step(rtdb, geom) ! Updates geom using ds 1613c 1614c We have now taken a step. Replace the approximate step taken 1615c by the exact step in case update of internals was not exact. 1616c 1617 call gsopt_compute_actual_step(geom) 1618c 1619c Energy and Gradient 1620c 1621 call mepgs_gra(rtdb, geom) 1622 if (.not. rtdb_get(rtdb,'task:energy', mt_dbl, 1, energy)) 1623 $ call errquit('gsopt: could not get energy',0, RTDB_ERR) 1624c 1625c Construct proyector 1626c 1627 call gsopt_cart_pmat(rtdb, geom) 1628 call dcopy(ncart, gx, 1, g, 1) ! g() set to gx() 1629c 1630c update the hessian 1631c 1632c 1633CCCCCCCCCCCCCCCCCCCCCC 1634CCCC MASS CCCC 1635CCCCCCCCCCCCCCCCCCCCCC 1636 if (mswg) then 1637 call mwcoord(ds, nvar, .true.) 1638 call mwgrad( g, nvar, .true.) 1639 call mwgrad( gp, nvar, .true.) 1640 end if 1641CCCCCCCCCCCCCCCCCCCCCC 1642CCCC MASS CCCC 1643CCCCCCCCCCCCCCCCCCCCCC 1644c 1645c 1646 call gsopt_hessian_update() 1647c 1648c 1649CCCCCCCCCCCCCCCCCCCCCC 1650CCCC MASS CCCC 1651CCCCCCCCCCCCCCCCCCCCCC 1652 if (mswg) then 1653 call mwcoord(ds, nvar, .false.) 1654 call mwgrad( g, nvar, .false.) 1655 call mwgrad( gp, nvar, .false.) 1656 end if 1657CCCCCCCCCCCCCCCCCCCCCC 1658CCCC MASS CCCC 1659CCCCCCCCCCCCCCCCCCCCCC 1660c 1661c 1662c Check time before next iteration 1663c 1664 required = int(1.2d0*(util_wallsec() - start)) + 1 1665 if (.not. util_test_time_remaining(rtdb,required)) goto 200 1666c 1667c 1668 enddo ! End of iterative loop 1669c istep = istep - 1 ! Since we fell out 1670 200 if (oprint.and.ga_nodeid().eq.0) write(6,201) 1671 201 format(/,1x,63('-')/,1x,'Failed to converge in maximum number', 1672 $ ' of steps or available time'/,1x,63('-')/) 1673 converged = .false. 1674c 1675 100 if (converged) then 1676 if (oprint.and.ga_nodeid().eq.0) write(6,101) 1677 101 format(/,6x,22('-'),/,6x,'Optimization converged',/, 1678 $ 6x,22('-'),/) 1679 endif 1680c 1681 if (ga_nodeid().eq.0 .and. util_print('finish',print_low)) then 1682c 1683c Print out final info and geometry 1684c 1685 call gsopt_print(geom, istep) 1686 if (.not. geom_print(geom)) call errquit 1687 $ ('hnd_opt_drv: geom_print?',0, GEOM_ERR) 1688c 1689c 1690 if (util_print('bonds',print_default)) then 1691 if (.not.geom_print_distances(geom)) call errquit( 1692 & 'hnd_opt_drv: geom_print_distances failed',911, 1693 & GEOM_ERR) 1694 endif 1695 if (util_print('angles',print_default)) then 1696 if (.not.geom_print_angles(geom)) call errquit( 1697 & 'hnd_opt_drv: geom_print_angles failed',911, 1698 & GEOM_ERR) 1699 endif 1700 endif 1701c 1702c Clean up and go home 1703c 1704 if (.not.geom_destroy(geom)) 1705 & call errquit('hnd_opt: geom_destroy?', 911, GEOM_ERR) 1706c 1707 gsopt=converged 1708 if (converged) then 1709 call ecce_print_module_exit('gsopt', 'ok') 1710 else 1711 call ecce_print_module_exit('gsopt', 'failed') 1712 endif 1713c 1714 call movecs_ecce_print_on() ! Restore MO printing 1715 call util_print_pop 1716c 1717 call ga_sync() 1718c 1719 end 1720c 1721 subroutine gsopt_hss_init(rtdb,geom) 1722 implicit none 1723#include "errquit.fh" 1724#include "mafdecls.fh" 1725#include "global.fh" 1726#include "geom.fh" 1727#include "rtdb.fh" 1728#include "nwc_const.fh" 1729#include "cgsopt.fh" 1730c 1731 integer rtdb, geom 1732 1733c 1734 double precision zero 1735 parameter (zero=0.0d+00) 1736 integer mxatom, mxcart, mxzmat, mxcoor 1737 parameter (mxatom=nw_max_atom) 1738 parameter (mxcart=3*mxatom) 1739 parameter (mxzmat=nw_max_zmat) 1740 parameter (mxcoor=nw_max_coor) 1741c 1742c These commons are used in the internal coordinate guess 1743c 1744 integer nuc 1745 COMMON/HND_MOLNUC/NUC(MXATOM) 1746 double precision c, zan 1747 integer natom 1748 common/hnd_molxyz/c(3,mxatom),zan(mxatom),natom 1749 integer nnzmat, nnzvar, nnvar 1750 common/hnd_zmtpar/nnzmat,nnzvar,nnvar 1751 double precision hscale, ascale, bscale, tscale, amat(3,3) 1752c 1753 integer l_hess, k_hess, l_zmat, k_zmat, l_izmat, k_izmat, i, j 1754 integer l_c, k_c, l_t, k_t, iat 1755 logical old_hessian 1756 character*16 atom_tags(mxatom) 1757c 1758 nnzmat = nzmat 1759 nnzvar = nzvar 1760 nnvar = nzvar 1761 if (.not. geom_ncent(geom,natom)) 1762 1 call errquit('hnd_opt: geom_ncent?',911, GEOM_ERR) 1763c 1764 if (.not. ma_push_get(mt_dbl, nvar**2, 'hessian', 1765 $ l_hess, k_hess)) call errquit 1766 $ ('gsopt_init_hess: failed allocating hessian',nvar**2, 1767 & MA_ERR) 1768c 1769 old_hessian=.false. 1770 if(inhess.ne.1) call gsopt_check_hess(nvar, old_hessian) 1771 if (oprint) write(6,*) 1772 if (old_hessian) then 1773 if (inhess.eq.2) then 1774 call gsopt_hess_cart_guess() 1775 if (oprint) write(6,*) 1776 $ ' Using Cartesian Hessian from previous frequency', 1777 $ ' calculation' 1778 else 1779 if (oprint) write(6,*) 'Using old Hessian from', 1780 $ ' previous optimization' 1781 end if 1782 goto 999 1783 else 1784 if (oprint) write(6,*) 'Using diagonal initial Hessian ' 1785 endif 1786c 1787c Cartesians are easy 1788c 1789 call dfill(nvar**2, 0.0d0, dbl_mb(k_hess), 1) 1790 call dfill(nvar, 0.5d0, dbl_mb(k_hess), nvar+1) 1791 if (isystype .ne. 0) then 1792 if (.not. geom_amatrix_get(geom, amat)) 1793 $ call errquit('geom_frac_to_cart: a', 0, GEOM_ERR) 1794 do iat = 1, nat_real 1795 do i = 1, 3 1796 dbl_mb(k_hess + (iat-1)*3 + (i-1) + 1797 $ ((iat-1)*3 + (i-1))*nat*3) = 0.5*amat(i,i)**2 1798 end do 1799 end do 1800 if (odebug) then 1801 write(6,*) ' The initial hessian ' 1802 call output(dbl_mb(k_hess),1,nvar,1,nvar,nvar,nvar,1) 1803 end if 1804 end if 1805c 1806c Artificially break degeneracies so that accidentally degenerate 1807c modes are split and therefore step restriction along modes 1808c is well defined 1809c 1810 do i = 1, nvar 1811 dbl_mb(k_hess+i-1 + (i-1)*nvar) = 1812 $ dbl_mb(k_hess+i-1 + (i-1)*nvar) + dble(i-1)*1e-7 1813 end do 1814c 1815 call geom_hnd_put_data('gsopt.hess', dbl_mb(k_hess), nvar*nvar) 1816c 1817c Apply constants, constraints and overall scaling 1818c 1819 999 if (.not. ma_push_get(mt_dbl, nvar**2, 'hessian', l_c, k_c)) 1820 $ call errquit('gsopt_init_hess: failed allocating hessian', 1821 $ nvar**2, MA_ERR) 1822 if (.not. ma_push_get(mt_dbl, nvar**2, 'hessian', l_t, k_t)) 1823 $ call errquit('gsopt_init_hess: failed allocating hessian', 1824 $ nvar**2, MA_ERR) 1825c 1826 call geom_hnd_get_data('gsopt.hess', dbl_mb(k_hess), nvar*nvar) 1827c 1828c Used to use c here ... now use p 1829c 1830 call geom_hnd_get_data('p',dbl_mb(k_c), nvar*nvar) 1831c 1832 if (odebug) then 1833 write(6,*) ' Initial Hessian before P' 1834 call output(dbl_mb(k_hess),1,nvar,1,nvar,nvar,nvar,1) 1835 endif 1836 call dgemm('n','n',nvar,nvar,nvar,1d0,dbl_mb(k_c),nvar, 1837 $ dbl_mb(k_hess),nvar,0d0,dbl_mb(k_t),nvar) 1838 call dgemm('n','t',nvar,nvar,nvar,1d0,dbl_mb(k_t),nvar, 1839 $ dbl_mb(k_c),nvar,0d0,dbl_mb(k_hess),nvar) 1840 if (odebug) then 1841 write(6,*) ' Initial Hessian after P' 1842 call output(dbl_mb(k_hess),1,nvar,1,nvar,nvar,nvar,1) 1843 endif 1844c 1845 if (nactive .ne. nat_real) then 1846c 1847c We are in cartesian coordinates and some have been frozen. 1848c Since there is no redundancy or coupling we just need 1849c to make sure that the initial Hessian does not couple 1850c frozen with unfrozen variables and we are OK. 1851c 1852 do iat = 1, nat 1853 if (.not. oactive(iat)) then 1854 do i = 1+(iat-1)*3, iat*3 1855 do j = 1, nvar 1856 dbl_mb(k_hess+j-1+(i-1)*nvar) = 0d0 1857 dbl_mb(k_hess+i-1+(j-1)*nvar) = 0d0 1858 enddo 1859 dbl_mb(k_hess+i-1+(i-1)*nvar) = 1d0 1860 enddo 1861 endif 1862 enddo 1863 endif 1864c 1865 if (.not. rtdb_get(rtdb,'gsopt:hscale',mt_dbl,1,hscale)) 1866 $ hscale = 1d0 1867 call dscal(nvar*nvar, hscale, dbl_mb(k_hess), 1) 1868 if (oprint .and. hscale.ne.1d0) then 1869 if (ga_nodeid().eq.0) write(6,78) hscale 1870 end if 1871 78 format(' Scaling initial hessian by ',f6.2) 1872c 1873 call geom_hnd_put_data('gsopt.hess',dbl_mb(k_hess), nvar*nvar) 1874c 1875 if (.not. ma_chop_stack(l_hess)) call errquit 1876 $ ('gsopt_init_hess ma corrupt',0, MA_ERR) 1877c 1878 end 1879 subroutine gsopt_check_hess(nvar, old_hessian) 1880 implicit none 1881#include "global.fh" 1882#include "tcgmsg.fh" 1883#include "mafdecls.fh" 1884c 1885 integer nvar 1886 logical old_hessian 1887 character*255 filename 1888c 1889 integer m 1890 double precision big,x 1891c 1892 big = 1d6 1893c 1894c Look at an existing hessian file and verify it 1895c 1896 call util_file_name('gsopt.hess',.false.,.false.,filename) 1897c 1898 if (ga_nodeid() .eq. 0) then 1899 open(32,file=filename,form='unformatted',status='old',err=10) 1900 read(32,err=11) m 1901 if (m.ne.nvar*nvar) goto 11 1902 close(32) 1903 old_hessian = .true. 1904 goto 20 1905c 1906 11 close(32) 1907 10 old_hessian = .false. 1908 endif 1909CJMC *** Correct for hessian reading from freq *** 1910 1911 20 call util_file_name('hess',.false.,.false.,filename) 1912c 1913 x = big 1914 if (ga_nodeid() .eq. 0) then 1915 open(32,file=filename,form='formatted',status='unknown', 1916 $ err=30,access='sequential') 1917 read(32,500,err=31,end=31) x 1918 if (x.eq.big) goto 31 1919 close(32) 1920 old_hessian = .true. 1921 goto 40 1922c 1923 31 close(32) 1924 30 old_hessian = .false. 1925 endif 1926CJMC *** Correct for hessian reading from freq *** 1927c 1928 40 call ga_brdcst(323, old_hessian, ma_sizeof(mt_log,1,mt_byte), 0) 1929c 1930 return 1931 500 format(f30.15) 1932 end 1933 subroutine gsopt_del_hess() 1934 implicit none 1935#include "util.fh" 1936c 1937c Delete the Hessian information restart file. 1938c 1939 character*255 opt_hess_fil 1940c 1941 call util_file_name('gsopt.hess', 1942 1 .false.,.false.,opt_hess_fil) 1943 call util_file_unlink(opt_hess_fil) 1944c 1945 if (util_print('information',print_low)) then 1946 write(6,*) 1947 write(6,*) ' Deleted TROPT restart files ' 1948 write(6,*) 1949 endif 1950c 1951 end 1952 subroutine gsopt_initialize(rtdb, geom) 1953 implicit none 1954#include "errquit.fh" 1955#include "nwc_const.fh" 1956#include "cmepgs.fh" 1957#include "cgsopt.fh" 1958#include "geom.fh" 1959#include "rtdb.fh" 1960#include "util.fh" 1961#include "global.fh" 1962#include "mafdecls.fh" 1963#include "inp.fh" 1964 integer rtdb 1965 integer geom ! [output] 1966c 1967c This routine initializes the common /cgsopt/ and 1968c also creates and returns the geometry handle 1969c 1970 integer i, j, num, ma_type, nactive_atoms, l_actlist 1971 integer ivar 1972 logical ignore 1973 character*80 title 1974 character*8 source, test 1975 character*32 theory 1976 logical gsopt_geom_cart_coords_get 1977c 1978 call util_print_push 1979 call util_print_rtdb_load(rtdb, 'gsopt') 1980 call ecce_print_module_entry('gsopt') 1981 oprint = util_print('information', print_low) 1982 $ .and. (ga_nodeid() .eq. 0) 1983 odebug = util_print('debug', print_debug) 1984 $ .and. (ga_nodeid() .eq. 0) 1985c 1986 if (rtdb_cget(rtdb,'title',1,title)) then 1987 if (oprint) then 1988 write(6,*) 1989 write(6,*) 1990 call util_print_centered(6, title, 40, .false.) 1991 write(6,*) 1992 write(6,*) 1993 endif 1994 endif 1995c 1996c ----- parameters for optimization gsopt ----- 1997c 1998 if (.not. rtdb_get(rtdb,'gsopt:modsad',mt_int,1,modsad)) 1999 $ modsad=0 2000cjmc 2001 trust = sqrt(ctrust2) 2002cjmc 2003 if (.not. rtdb_cget(rtdb,'mepgs:xyz',1,xyz)) 2004 $ xyz = ' ' 2005 if (.not. rtdb_get(rtdb,'gsopt:eprec',mt_dbl,1,eprec)) then 2006 if (.not. rtdb_cget(rtdb,'task:theory',1,theory)) 2007 $ theory = ' ' 2008 if (inp_compare(.false.,theory,'dft')) then 2009 eprec = 5e-6 2010 else 2011 eprec = 1e-7 2012 endif 2013 endif 2014 if (.not. rtdb_get(rtdb,'gsopt:gmax_tol',mt_dbl,1,gmax_tol)) 2015 $ gmax_tol = 0.00045d0 2016 if (.not. rtdb_get(rtdb,'gsopt:grms_tol',mt_dbl,1,grms_tol)) 2017 $ grms_tol = 0.0003d0 2018 if (.not. rtdb_get(rtdb,'gsopt:xmax_tol',mt_dbl,1,xmax_tol)) 2019 $ xmax_tol = 0.0018d0 2020 if (.not. rtdb_get(rtdb,'gsopt:xrms_tol',mt_dbl,1,xrms_tol)) 2021 $ xrms_tol = 0.0012d0 2022 if (.not. rtdb_get(rtdb,'gsopt:nptopt',mt_int,1,nptopt)) 2023 $ nptopt=50 2024 if (.not. rtdb_get(rtdb,'gsopt:inhess',mt_int,1,inhess)) 2025 $ inhess=0 2026 if (.not. rtdb_get(rtdb,'gsopt:linopt',mt_int,1,linopt)) 2027 $ linopt=1 2028 if (.not. rtdb_get(rtdb,'gsopt:moddir',mt_int,1,moddir)) 2029 $ moddir=0 2030 if (.not. rtdb_get(rtdb,'gsopt:modsad',mt_int,1,modsad)) 2031 $ modsad=0 2032 if (.not. rtdb_get(rtdb,'gsopt:sadstp',mt_dbl,1,sadstp)) 2033 $ sadstp=0.1d0 2034 if (.not. rtdb_get(rtdb,'gsopt:oqstep',mt_log,1,oqstep)) 2035 $ oqstep = .true. 2036 if (.not. rtdb_get(rtdb,'gsopt:modupd',mt_int,1,modupd)) then 2037 if (modsad .eq. 0) then 2038 modupd = 1 ! BFGS update for minimization 2039 else 2040 modupd = 2 ! PSB update for saddle point 2041 endif 2042 endif 2043 if (.not. rtdb_get(rtdb,'gsopt:ocheckgrad',mt_log,1,ocheckgrad)) 2044 $ ocheckgrad = .false. 2045 2046 if (.not. rtdb_get(rtdb,'includestress',mt_log,1,ostress)) then 2047 ostress = .false. 2048 end if 2049 if (.not. rtdb_get(rtdb,'includelattice',mt_log,1,ostress2)) then 2050 ostress2 = .false. 2051 end if 2052 if ((.not.ostress).and.(ostress2)) ostress2 = .false. 2053 if ((ostress) .and.(ostress2)) ostress = .false. 2054 2055c 2056c Force sensible options 2057c 2058 if (modsad .eq. 0) then 2059 modupd = 1 ! BFGS update for minimization 2060 else 2061 linopt = 0 ! No line search for saddle 2062 endif 2063c 2064c Save a copy of the initial geometry so we can analyze what 2065c happened during the optimization 2066c 2067 if (ga_nodeid() .eq. 0) then 2068 ignore = rtdb_parallel(.false.) 2069 if (.not. geom_create(geom, 'geometry')) 2070 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 2071 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 2072 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 2073 if (.not. geom_rtdb_store(rtdb, geom, 'gsoptinitial')) 2074 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 2075 if (.not. geom_destroy(geom)) 2076 $ call errquit('gsopt: geom_destroy?',0, GEOM_ERR) 2077 ignore = rtdb_parallel(.true.) 2078 endif 2079 call ga_sync() 2080c 2081c Load the geometry info 2082c 2083 if (.not. geom_create(geom, 'geometry')) 2084 & call errquit('hnd_opt: geom_create?', 911, GEOM_ERR) 2085 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 2086 & call errquit('hnd_opt: no geometry ', 911, RTDB_ERR) 2087 if (.not. geom_ncent(geom,nat)) 2088 $ call errquit('hnd_opt: natoms?',nat, GEOM_ERR) 2089 call grad_active_atoms(rtdb, nat, oactive, nactive) 2090 if (.not. geom_systype_get(geom, isystype)) 2091 $ call errquit('gsopt: systype?',0, GEOM_ERR) 2092c 2093 if (oprint.and.ga_nodeid().eq.0) then 2094 write(6,1) gmax_tol, grms_tol, xmax_tol, xrms_tol, 2095 $ eprec, 2096 $ nptopt, inhess, modupd 2097 1 format( 2098 $ ' maximum gradient threshold (gmax) = ', f10.6,/, 2099 $ ' rms gradient threshold (grms) = ', f10.6,/, 2100 $ ' maximum cartesian step threshold (xmax) = ', f10.6,/, 2101 $ ' rms cartesian step threshold (xrms) = ', f10.6,/, 2102 $ ' energy precision (eprec) = ', 1p,d9.1, 2103 $ 0p,/, 2104 $ ' maximum number of steps (nptopt) = ', i4,/, 2105 $ ' initial hessian option (inhess) = ', i4,/, 2106 $ ' hessian update option (modupd) = ', i4,/) 2107 if (modsad .eq. 0) then 2108 write(6,9994) 2109 9994 format(/,10x,19('-'), 2110 1 /,10x,'Energy Minimization', 2111 2 /,10x,19('-'),/) 2112 else 2113 write(6,9995) 2114 9995 format(/,10x,23('-'), 2115 1 /,10x,'Transition State Search', 2116 2 /,10x,23('-'),/) 2117 endif 2118 if (ostress) then 2119 write(6,*) ' INCLUDING STRESS !!!!!!!!!!!!!!!!' 2120 if (isystype.eq.0) call errquit('NOT A PERIODIC SYSTEM',0, 2121 & GEOM_ERR) 2122 endif 2123 if (ostress2) then 2124 write(6,*) ' INCLUDING LATTICE GRADIENTS !!!!!' 2125 if (isystype.eq.0) call errquit('NOT A PERIODIC SYSTEM',0, 2126 & GEOM_ERR) 2127 endif 2128 2129 call util_flush(6) 2130 endif 2131c 2132c Nvar is the no. of variables in the optimization 2133c 2134c If we are optimizing the unit cell parameters then we pretend 2135c there there are 3 more atoms which will parameterize the 2136c unit cell. 2137c 2138 nat_real = nat 2139 if (ostress) nat = nat + 3 2140 if (ostress2) nat = nat + 2 2141 ncart = 3*nat 2142 nvar = ncart 2143c 2144c 2145 call gsopt_cart_pmat(rtdb, geom) 2146c 2147 energy = 0d0 2148 energyp= 0d0 2149 alpha = 1d0 2150 gmax = 0d0 2151 grms = 0d0 2152 smax = 0d0 2153 srms = 0d0 2154 xmax = 0d0 2155 xrms = 0d0 2156 call dfill(max_nvar, 0d0, ds, 1) 2157 call dfill(max_nvar, 0d0,dsp, 1) 2158 call dfill(max_nvar, 0d0, gx, 1) 2159 call dfill(max_nvar, 0d0, gq, 1) 2160 call dfill(max_nvar, 0d0, g, 1) 2161 call dfill(max_nvar, 0d0, gp, 1) 2162 call dfill(max_nvar, 0d0, radius, 1) 2163c 2164 if (.not. gsopt_geom_cart_coords_get(geom, sp)) 2165 $ call errquit('gsopt: geom?',0, GEOM_ERR) 2166 call dcopy(nvar, sp, 1, radius, 1) 2167 call daxpy(nvar, -1.0d0, center, 1, radius, 1) 2168c 2169 end 2170 subroutine gsopt_hessian_update() 2171 implicit none 2172#include "errquit.fh" 2173#include "nwc_const.fh" 2174#include "cmepgs.fh" 2175#include "cgsopt.fh" 2176#include "util.fh" 2177#include "mafdecls.fh" 2178c 2179c Update the current Hessian in the optimization variables using 2180c . gp() - the gradient at the previous point 2181c . g() - the gradient at the current point 2182c . ds() - the previous search direction 2183c . alpha - the step in the previous search direction 2184c 2185c Only the Hessian is modified. 2186c 2187 double precision hds(max_nvar) 2188 double precision dsds, dshds, dsdg 2189 integer l_hess, k_hess, i, j 2190 integer ind 2191 ind(i,j) = k_hess + i + (j-1)*nvar - 1 2192c 2193 if (alpha .eq. 0d0) call errquit 2194 $ ('gsopt_hessian_update: zero step?',0, GEOM_ERR) 2195 call dscal(nvar, alpha, ds, 1) 2196c 2197 if (.not. ma_push_get(mt_dbl, nvar**2, 'hess', 2198 $ l_hess, k_hess)) call errquit 2199 $ ('gsopt_hessian_update: memory for hessian',nvar**2, 2200 & GEOM_ERR) 2201 call geom_hnd_get_data('gsopt.hess',dbl_mb(k_hess), nvar**2) 2202c 2203c Form bits and pieces that are needed 2204c 2205 call dgemv('n',nvar,nvar,1d0,dbl_mb(k_hess),nvar, 2206 $ ds,1,0d0,hds,1) 2207c 2208 dshds = ddot(nvar, ds, 1, hds, 1) 2209 dsds = ddot(nvar, ds, 1, ds, 1) 2210 dsdg = 0d0 2211 do i = 1, nvar 2212 dsdg = dsdg + ds(i)*(g(i) - gp(i)) 2213 enddo 2214c 2215 if(modupd.le.1) then 2216c 2217c ----- -bfgs- update ----- 2218c 2219 if(dsdg.gt.1d-14) then 2220 do i=1,nvar 2221 do j=1,nvar 2222 dbl_mb(ind(i,j))=dbl_mb(ind(i,j)) 2223 $ + (g(i)-gp(i))*(g(j)-gp(j))/dsdg 2224 1 - hds(i)* hds(j)/dshds 2225 enddo 2226 enddo 2227 endif 2228 else 2229c 2230c ----- -psb- update ----- 2231c 2232 if (abs(dsdg).gt.1d-8) then 2233 do i=1,nvar 2234 do j=1,nvar 2235 dbl_mb(ind(i,j))=dbl_mb(ind(i,j)) 2236 $ + ((g(i)-gp(i))-hds(i))*ds(j)/dsds 2237 1 + ((g(j)-gp(j))-hds(j))*ds(i)/dsds 2238 2 - ds(i)*ds(j)*(dsdg-dshds)/(dsds*dsds) 2239 enddo 2240 enddo 2241 endif 2242 endif 2243c 2244 call geom_hnd_put_data('gsopt.hess',dbl_mb(k_hess), nvar**2) 2245 if (.not. ma_pop_stack(l_hess)) call errquit 2246 $ ('gsopt_hessian_update: ma?',0, MA_ERR) 2247c 2248c 2249 end 2250c 2251 subroutine gsopt_take_step(rtdb, geom) 2252 implicit none 2253#include "errquit.fh" 2254#include "nwc_const.fh" 2255#include "cmepgs.fh" 2256#include "cgsopt.fh" 2257#include "geom.fh" 2258#include "rtdb.fh" 2259#include "mafdecls.fh" 2260 integer rtdb, geom 2261c 2262c Update the geometry in geom and in the database 2263c 'geometry' by taking the step 2264c alpha*ds() in the optimization variables 2265c 2266c The geom is modified, and xmax/xrms are computed from the 2267c first-order step. 2268c 2269 double precision xold(max_cart), xnew(max_cart), err 2270 double precision aaa(3,3) 2271 integer i, l_bi, k_bi 2272 logical gsopt_geom_cart_coords_get 2273 logical gsopt_geom_cart_coords_set 2274c 2275c Get original coordinates 2276c 2277c FRACTIONAL? 2278 if (.not. gsopt_geom_cart_coords_get(geom, xold)) 2279 $ call errquit('gsopt_energy_step: coordinates?',geom, 2280 & GEOM_ERR) 2281c 2282c Take the step 2283c 2284 call dcopy(ncart, ds, 1, xnew, 1) 2285 call dscal(nvar, alpha, xnew, 1) 2286 call sym_grad_symmetrize(geom, xnew) 2287 call daxpy(ncart, 1.0d0, xold, 1, xnew, 1) 2288c FRACTIONAL? 2289 if (.not. gsopt_geom_cart_coords_set(geom, xnew)) 2290 $ call errquit('gsopt_energy_step: coordinates?',geom, 2291 & GEOM_ERR) 2292c 2293c Must ensure the geometry has the required symmetry even after 2294c enforcing it on the step. Should use 2295c an error criterion consistent with the step size. 2296c 2297 call sym_geom_project(geom, trust) 2298c 2299 if (.not. geom_rtdb_store(rtdb, geom, 'geometry')) 2300 $ call errquit('gsopt_energy_step: grs?',geom, RTDB_ERR) 2301c 2302c Compute the maximum and RMS cartesian displacements 2303c 2304c 2305CCCCCCCCCCCCCCCCCCCCCC 2306CCCC MASS CCCC 2307CCCCCCCCCCCCCCCCCCCCCC 2308 if (mswg) then 2309 call mwcoord(xold, nvar, .true.) 2310 call mwcoord(xnew, nvar, .true.) 2311 end if 2312CCCCCCCCCCCCCCCCCCCCCC 2313CCCC MASS CCCC 2314CCCCCCCCCCCCCCCCCCCCCC 2315c 2316c 2317 xmax = 0d0 2318 xrms = 0d0 2319 do i = 1, ncart 2320 xmax = max(xmax, abs(xold(i)-xnew(i))) 2321 xrms = xrms + (xold(i)-xnew(i))**2 2322 enddo 2323 xrms = sqrt(xrms/dble(ncart)) 2324c 2325c 2326CCCCCCCCCCCCCCCCCCCCCC 2327CCCC MASS CCCC 2328CCCCCCCCCCCCCCCCCCCCCC 2329 if (mswg) then 2330 call mwcoord(xold, nvar, .false.) 2331 call mwcoord(xnew, nvar, .false.) 2332 end if 2333CCCCCCCCCCCCCCCCCCCCCC 2334CCCC MASS CCCC 2335CCCCCCCCCCCCCCCCCCCCCC 2336c 2337c 2338 end 2339 subroutine gsopt_print(geom, istep) 2340 implicit none 2341#include "errquit.fh" 2342#include "nwc_const.fh" 2343#include "cgsopt.fh" 2344#include "global.fh" 2345#include "util.fh" 2346#include "inp.fh" 2347 2348 integer geom, istep 2349c 2350c Print out stuff 2351c 2352 integer i 2353 double precision de 2354 character*9 cvg1, cvg2, cvg3, cvg4 2355 character*1 mark 2356c 2357 de = 0d0 2358 if (istep .gt. 1) de = energy-energyp 2359 cvg1 = ' ' 2360 cvg2 = ' ' 2361 cvg3 = ' ' 2362 cvg4 = ' ' 2363 if (gmax .lt. gmax_tol) cvg1 = ' ok ' 2364 if (grms .lt. grms_tol) cvg2 = ' ok ' 2365 if (xrms .lt. xrms_tol) cvg3 = ' ok ' 2366 if (xmax .lt. xmax_tol) cvg4 = ' ok ' 2367c 2368 if (oprint) then 2369 mark = '@' 2370 if (istep .gt. 1) mark = ' ' 2371 if (ga_nodeid().eq.0) write(6,1) mark, mark 2372 mark = '@' 2373 if (ga_nodeid().eq.0) write(6,2) mark, istep-1, energy, de, 2374 $ gmax, grms, xrms, xmax, util_wallsec(), 2375 $ cvg1, cvg2, cvg3, cvg4 2376 1 format( 2377 $ /,a1,' Step Energy Delta E Gmax', 2378 $ ' Grms Xrms Xmax Walltime', 2379 $ /,a1,' ---- ---------------- -------- --------', 2380 $ ' -------- -------- -------- --------') 2381 2 format( 2382 $ a1,i5,f17.8,1p,d9.1,0p,4f9.5,f9.1,/, 2383 $ 1x,5x,17x,9x,4a9,/) 2384 endif 2385c 2386 end 2387 logical function gsopt_converged(optcyc) 2388 implicit none 2389#include "nwc_const.fh" 2390#include "cgsopt.fh" 2391c 2392 integer optcyc 2393 double precision de 2394c 2395c Return true if we have converged 2396c 2397c Nothing is modified. Assumes gsopt_compute_info() 2398c has been called. 2399c gmax_tol, ! [user] tolerance for max internal gradient 2400c grms_tol, ! [user] tolerance for rms internal gradient 2401c xrms_tol, ! [user] tolerance for rms cartesian step 2402c xmax_tol, ! [user] tolerance for max cartesian step 2403c 2404 de = abs(energy-energyp) 2405 gsopt_converged = 2406 $ ((gmax .lt. gmax_tol) .and. (grms .lt. grms_tol) .and. 2407 $ (xrms .lt. xrms_tol) .and. (xmax .lt. xmax_tol)) 2408 $ .or. 2409 $ ((gmax.lt.0.01d0*gmax_tol) .and. (grms.lt.0.01d0*grms_tol)) 2410CJMC 2411c 2412c *** Convergence over energy criterion 2413c 2414 $ .or. 2415 $ ((gmax .lt. gmax_tol) .and. (grms .lt. grms_tol).and. 2416 $ (de .lt. eprec)) 2417 $ .or. 2418c 2419c *** Initial structure already a critical point 2420c 2421 $ ((gmax .lt. gmax_tol) .and. (grms .lt. grms_tol).and. 2422 $ (optcyc .eq. 1)) 2423CJMC 2424 if ((gmax .lt. gmax_tol) .and. (grms .lt. grms_tol).and. 2425 $ (de .lt. eprec)) then 2426 write(6,*) "Convergence reached over gradient and energy" 2427 write(6,*) "Energy change = " , de 2428 write(6,*) "Energy precision = " , eprec 2429 end if 2430 2431c 2432 end 2433 subroutine gsopt_cart_pmat(rtdb, geom) 2434 implicit none 2435#include "errquit.fh" 2436#include "rtdb.fh" 2437#include "nwc_const.fh" 2438#include "cgsopt.fh" 2439#include "geom.fh" 2440#include "mafdecls.fh" 2441#include "util.fh" 2442 integer rtdb, geom 2443c 2444c Compute the cartesian equivalent of the P = G.G^-1 matrix 2445c which projects to and from the linearly independent 2446c set of coordinates. In the cartesian case P is the complement 2447c of the projector onto the rotations and translations 2448c For ease of use we also write out a unit Binv matrix. 2449c 2450c Only the P/Binv matrices are generated. Nothing is modified. 2451c 2452c Minor little catch is that if some atoms are being frozen 2453c we are no longer invariant to translations or rotations. 2454c 2455c RTDB is used to look for frozen atoms. 2456c 2457 double precision centroid(3), x, y, z, xx, yy, zz, fx 2458 double precision coords(3,max_cent) 2459 double precision work(max_cart,6) 2460 integer i, j, k, l_pmat, k_pmat, i3, ma_type, nelem 2461 character*26 date 2462 integer ind 2463 logical task_qmmm 2464 logical gsopt_geom_cart_coords_get 2465 ind(i,j) = k_pmat + i-1 + (j-1)*ncart 2466c 2467c FRACTIONAL? 2468 if (.not. gsopt_geom_cart_coords_get(geom, coords)) 2469 $ call errquit('gsopt_cart_pmat: geom?',geom, GEOM_ERR) 2470c 2471c Construct normalized vectors in work in the direction 2472c of the rotations and translations. 2473c 2474 call dfill(3, 0.0d0, centroid, 1) 2475 do i = 1, nat 2476 do k = 1, 3 2477 centroid(k) = centroid(k) + coords(k,i)/nat 2478 enddo 2479 enddo 2480c 2481 do k = 1, 3 ! x, y, z translations 2482 call dfill(ncart, 0.0d0, work(1,k), 1) 2483 call dfill(nat, sqrt(1.0d0/nat), work(k,k), 3) 2484 enddo 2485 do k = 4, 6 ! x, y, z rotations 2486 do i = 1, nat 2487 x = coords(1,i) - centroid(1) 2488 y = coords(2,i) - centroid(2) 2489 z = coords(3,i) - centroid(3) 2490 if (k .eq. 4) then 2491 xx = 0.0d0 2492 yy = -z 2493 zz = y 2494 else if (k .eq. 5) then 2495 xx = z 2496 yy = 0.0d0 2497 zz = -x 2498 else if (k .eq. 6) then 2499 xx = -y 2500 yy = x 2501 zz = 0.0d0 2502 endif 2503 i3 = (i-1)*3 2504 work(i3+1,k) = xx 2505 work(i3+2,k) = yy 2506 work(i3+3,k) = zz 2507 enddo 2508 do j = 1, k-1 2509 fx = ddot(ncart, work(1,j), 1, work(1,k), 1) 2510 call daxpy(ncart, -fx, work(1,j), 1, work(1,k), 1) 2511 enddo 2512 fx = sqrt(ddot(ncart, work(1,k), 1, work(1,k), 1)) 2513 if (fx . gt. 1d-6) then 2514 call dscal(ncart, 1.0d0/fx, work(1,k), 1) 2515 else 2516 call dfill(ncart, 0.0d0, work(1,k), 1) 2517 endif 2518 enddo 2519c 2520c The project is then 1 - V.VT where V is in work 2521c 2522 if (.not. ma_push_get(mt_dbl, ncart**2, 'pmat', 2523 $ l_pmat, k_pmat)) call errquit 2524 $ ('gsopt_cart_pmat: memory for pmat',ncart**2, GEOM_ERR) 2525c 2526c Form unit matrix 2527c 2528 call dfill(ncart**2, 0d0, dbl_mb(k_pmat), 1) 2529 call dfill(ncart, 1d0, dbl_mb(k_pmat), ncart+1) 2530c 2531c Store dummy unit matrix for B, Binv ... the cartesian 2532c gradient should already be invariant to rotations and translations. 2533c Also store dummy unit matrix for cmat (constraints) 2534c 2535 call geom_hnd_put_data('b', dbl_mb(k_pmat), ncart**2) 2536 call geom_hnd_put_data('b^-1', dbl_mb(k_pmat), ncart**2) 2537 call geom_hnd_put_data('c', dbl_mb(k_pmat), ncart**2) 2538c 2539 if (.not.rtdb_get(rtdb,'task:QMMM',mt_log,1,task_qmmm)) 2540 & task_qmmm = .false. 2541 2542 if ( rtdb_get_info(rtdb, 'geometry:actlist', ma_type, 2543 $ nelem, date) .or. 2544 $ rtdb_get_info(rtdb, 'geometry:inactlist', ma_type, 2545 $ nelem, date) .or. 2546 $ isystype .ne. 0 .or. 2547 $ task_qmmm .or. 2548 $ geom_extbq_on() ) then 2549c 2550c Some atoms are frozen or we have a periodic system so don't have 2551c invariance ... also store unit matrix for P. 2552c or we have QMMM calculation here 2553c 2554 call geom_hnd_put_data('p', dbl_mb(k_pmat), ncart**2) 2555c 2556 else 2557c 2558c Finish P 2559c 2560 do i = 1, ncart 2561 do j = 1, ncart 2562 do k = 1, 6 2563 dbl_mb(ind(j,i)) = dbl_mb(ind(j,i)) - 2564 $ work(j,k)*work(i,k) 2565 enddo 2566 enddo 2567 enddo 2568c 2569 call geom_hnd_put_data('p', dbl_mb(k_pmat), ncart**2) 2570c 2571 if (odebug) then 2572 write(6,*) ' Cartesian P matrix' 2573 call output(dbl_mb(k_pmat),1,ncart,1,ncart,ncart,ncart,1) 2574 endif 2575 endif 2576c 2577 if (.not. ma_chop_stack(l_pmat)) call errquit 2578 $ ('gsopt_cart_bmat: ma?',0, MA_ERR) 2579c 2580 end 2581 subroutine gsopt_project_hess_grad(hess, pg) 2582 implicit none 2583#include "errquit.fh" 2584#include "nwc_const.fh" 2585#include "cgsopt.fh" 2586#include "mafdecls.fh" 2587 double precision 2588 $ hess(nvar,nvar), ! returns projected & shifted Hessian 2589 $ pg(nvar) ! returns projected gradient 2590c 2591c Project and shift the Hessian and gradient following Peng et al. 2592c 2593c Nothing else is changed. 2594c 2595 integer l_pmat, k_pmat, l_work, k_work, i 2596 double precision big 2597c 2598 if (.not. ma_push_get(mt_dbl, nvar**2, 'work', 2599 $ l_work, k_work)) call errquit 2600 $ ('gsopt_proj_h_g: memory for pmat',nvar**2, MA_ERR) 2601 if (.not. ma_push_get(mt_dbl, nvar**2, 'pmat', 2602 $ l_pmat, k_pmat)) call errquit 2603 $ ('gsopt_proj_h_g: memory for work',nvar**2, MA_ERR) 2604c 2605 call geom_hnd_get_data('gsopt.hess',hess, nvar**2) 2606 if (odebug) then 2607 write(6,*) ' Hessian before projection' 2608 call output(hess, 1, nvar, 1, nvar, nvar, nvar, 1) 2609 write(6,*) ' Gradient before projection' 2610 call doutput(g, 1, nvar, 1, 1, nvar, 1, 1) 2611 endif 2612 call geom_hnd_get_data('p',dbl_mb(k_pmat), nvar**2) 2613 if (.not. ma_verify_allocator_stuff()) 2614 $ call errquit('freddy',0, MA_ERR) 2615c 2616c PG 2617c 2618 call dgemv('n',nvar, nvar, 1d0, dbl_mb(k_pmat), nvar, 2619 $ g, 1, 0d0, pg, 1) 2620 if (odebug) then 2621 write(6,*) ' Gradient after projection' 2622 call doutput(g, 1, nvar, 1, 1, nvar, 1, 1) 2623 endif 2624c 2625c PHP + 1000*(1-P) 2626c 2627 call dgemm('n', 'n', nvar, nvar, nvar, 1d0, dbl_mb(k_pmat), nvar, 2628 $ hess, nvar, 0d0, dbl_mb(k_work), nvar) 2629 call dgemm('n', 'n', nvar, nvar, nvar, 1d0, dbl_mb(k_work), nvar, 2630 $ dbl_mb(k_pmat), nvar, 0d0, hess, nvar) 2631 if (odebug) then 2632 write(6,*) ' Hessian after projection before shift' 2633 call output(hess, 1, nvar, 1, nvar, nvar, nvar, 1) 2634 endif 2635c 2636CJMC big = 1000d0 2637CJMC call daxpy(nvar*nvar, -big, dbl_mb(k_pmat), 1, hess, 1) 2638CJMC do i = 1, nvar 2639CJMC hess(i,i) = hess(i,i) + big 2640CJMC enddo 2641 if (odebug) then 2642 write(6,*) ' Hessian after projection ' 2643 call output(hess, 1, nvar, 1, nvar, nvar, nvar, 1) 2644 endif 2645c 2646 if (.not. ma_chop_stack(l_work)) call errquit 2647 $ ('gsopt_p_h_g:ma?',0, MA_ERR) 2648c 2649 end 2650 subroutine gsopt_compute_info() 2651 implicit none 2652#include "nwc_const.fh" 2653#include "cmepgs.fh" 2654#include "cgsopt.fh" 2655#include "util.fh" 2656 double precision desphere(max_nvar) 2657 double precision zeta(max_nvar) 2658 double precision norm 2659c 2660c Compute stuff used for printing and convergence tests 2661c 2662c gmax = maxmimum gradient element in optimization variables 2663c grms = rms grad 2664c smax = maximum step in opt. var 2665c srms = rms step 2666c 2667c xrms and xmax are computed by gsopt_take_step from the 2668c first-order step. 2669c 2670 integer i 2671c 2672 grms = 0d0 2673 srms = 0d0 2674 gmax = 0d0 2675 smax = 0d0 2676c 2677 call dfill(nvar, 0d0, zeta, 1) 2678 call dcopy(nvar, radius, 1, zeta, 1) 2679 call daxpy(nvar, 1.0d00, ds, 1, zeta, 1) 2680 norm = ddot(nvar, zeta, 1, zeta, 1) 2681 norm = ddot(nvar, g, 1, zeta, 1)/norm 2682 call dfill(nvar, 0d0, desphere, 1) 2683 call dcopy(nvar, g, 1, desphere, 1) 2684 call daxpy(nvar, -norm, zeta, 1, desphere, 1) 2685c 2686 do i = 1, nvar 2687 grms = grms + desphere(i)*desphere(i) 2688 srms = srms + ds(i)*ds(i) 2689 gmax = max(gmax, abs(desphere(i))) 2690 smax = max(smax, abs(ds(i))) 2691 enddo 2692 grms = sqrt(grms/dble(nvar)) 2693 srms = sqrt(srms/dble(nvar)) 2694c 2695 end 2696 subroutine gsopt_hess_cart_guess() 2697 implicit none 2698#include "errquit.fh" 2699#include "mafdecls.fh" 2700#include "global.fh" 2701#include "nwc_const.fh" 2702#include "cgsopt.fh" 2703#include "inp.fh" 2704c 2705c Read in cartesian Hessian and transform it as necessary 2706c to internal coordinates (neglecting the component due to 2707c the derivative) and writing the result to the hessian file. 2708c 2709c Reads file in vib_vib format using vib_vib filename default 2710c Note the default filename is set in task_freq 2711c filenames must be made identical. 2712c 2713c Format of vib file is ascii lower triangular elements only. 2714c 2715 integer h_unit 2716 parameter (h_unit=47) 2717 character*255 fname 2718 double precision x 2719 integer i,j 2720 integer l_bi, k_bi, l_hc, k_hc, l_hq, k_hq 2721c 2722 if (.not. ma_push_get(mt_dbl, ncart*nvar, 'binv', 2723 $ l_bi, k_bi)) call errquit 2724 $ ('gsopt_hess_cart_guess: ma?', ncart*nvar, MA_ERR) 2725c 2726 if (.not. ma_push_get(mt_dbl, max(ncart**2,nvar**2), 'hcart', 2727 $ l_hc, k_hc)) call errquit 2728 $ ('gsopt_hess_cart_guess: ma?', ncart**2, MA_ERR) 2729c 2730 if (.not. ma_push_get(mt_dbl, max(ncart**2,nvar**2), 'hcart2', 2731 $ l_hq, k_hq)) call errquit 2732 $ ('gsopt_hess_cart_guess: ma?', nvar**2, MA_ERR) 2733c 2734 if (ga_nodeid().eq.0) then 2735 call util_file_name('hess',.false.,.false.,fname) 2736 open(unit=h_unit,file=fname,form='formatted',status='unknown', 2737 $ err=99990,access='sequential') 2738 rewind h_unit 2739 do i = 1,ncart 2740 do j = 1,i 2741 read(h_unit,10000,err=99992,end=99992) x 2742 dbl_mb(k_hc+(i-1)*ncart+(j-1)) = x 2743 dbl_mb(k_hc+(j-1)*ncart+(i-1)) = x 2744 enddo 2745 enddo 2746 close(unit=h_unit,status='keep') 2747 endif 2748 call ga_brdcst(1,dbl_mb(k_hc),8*ncart**2,0) 2749c 2750 call geom_hnd_get_data('b^-1', dbl_mb(k_bi), nvar*ncart) 2751 call dgemm('n', 'n', ncart, nvar, ncart, 1d0, dbl_mb(k_hc), ncart, 2752 $ dbl_mb(k_bi), ncart, 0d0, dbl_mb(k_hq), ncart) 2753 call dgemm('t', 'n', nvar, nvar, ncart, 1d0, dbl_mb(k_bi), ncart, 2754 $ dbl_mb(k_hq), ncart, 0d0, dbl_mb(k_hc), nvar) 2755c 2756 do i = 1,nvar 2757 do j = 1,i 2758 x = (dbl_mb(k_hc+(i-1)*nvar+(j-1)) + 2759 $ dbl_mb(k_hc+(j-1)*nvar+(i-1))) * 0.5d0 2760 dbl_mb(k_hc+(i-1)*nvar+(j-1)) = x 2761 dbl_mb(k_hc+(j-1)*nvar+(i-1)) = x 2762 enddo 2763 enddo 2764c 2765 call geom_hnd_put_data('gsopt.hess',dbl_mb(k_hc), nvar**2) 2766c 2767 if (.not. ma_chop_stack(l_bi)) 2768 $ call errquit('gsopt_hess_cart_guess: ma corrupt?',0, MA_ERR) 2769c 2770 return 277110000 format(f30.15) 277299990 write(6,*)' could not open <',fname(1:inp_strlen(fname)), 2773 $ '> as unknown file' 2774 call errquit('gsopt_hess_cart: fatal error', 911, GEOM_ERR) 277599991 write(6,*)' could not open <',fname(1:inp_strlen(fname)), 2776 $ '> as new file' 2777 call errquit('gsopt_hess_cart: fatal error', 911, GEOM_ERR) 277899992 write(6,*)' error in reading <',fname(1:inp_strlen(fname)), 2779 $ '> as hessian file' 2780 call errquit('gsopt_hess_cart: fatal error', 911, GEOM_ERR) 2781 end 2782 subroutine gsopt_symmetrize_step(geom) 2783 implicit none 2784#include "errquit.fh" 2785#include "nwc_const.fh" 2786#include "cgsopt.fh" 2787#include "mafdecls.fh" 2788 integer geom 2789c 2790c Force symmetry upon the current search direction by projecting 2791c ds() onto symmetric component in cartesians 2792c 2793c Updates ds(). 2794c 2795 double precision dx(max_nvar) 2796 integer k_bi, l_bi 2797c 2798c 1) dx = dq*B-1 2799c 2) symmetrize dx 2800c 3) dq = dx*B 2801c 2802 if (ostress) return ! Yikes! 2803 if (ostress2) return ! Yikes! 2804c 2805 if (.not. ma_push_get(mt_dbl, ncart*nvar,'binv',l_bi, k_bi)) 2806 $ call errquit('gsopt_sym: memory for binv', ncart*nvar, 2807 & MA_ERR) 2808 call geom_hnd_get_data('b^-1', dbl_mb(k_bi), ncart*nvar) 2809 if (odebug) then 2810 write(6,*) ' Symmetrize step - initial q ' 2811 call doutput(ds, 1, nvar, 1, 1, nvar, 1, 1) 2812 endif 2813 call dgemv('n', ncart, nvar, 1d0, dbl_mb(k_bi), ncart, 2814 $ ds, 1, 0.0d0, dx, 1) 2815 call sym_grad_symmetrize(geom, dx) 2816 call geom_hnd_get_data('b', dbl_mb(k_bi), ncart*nvar) 2817 call dgemv('t', ncart, nvar, 1d0, dbl_mb(k_bi), ncart, 2818 $ dx, 1, 0.0d0, ds, 1) 2819 if (odebug) then 2820 write(6,*) ' Symmetrize step - final q ' 2821 call doutput(ds, 1, nvar, 1, 1, nvar, 1, 1) 2822 endif 2823c 2824 if (.not. ma_pop_stack(l_bi)) call errquit('gsopt_sym:ma',0, 2825 & MA_ERR) 2826c 2827 end 2828 subroutine gsopt_compute_actual_step(geom) 2829 implicit none 2830#include "errquit.fh" 2831#include "nwc_const.fh" 2832#include "cmepgs.fh" 2833#include "cgsopt.fh" 2834#include "geom.fh" 2835 integer geom 2836c 2837c We have now taken a step. Since the non-linear transformations 2838c involved in taking a step may not have been done exactly replace 2839c ds() with the actual step taken so that the Hessian may be precisely 2840c updated. Updates ds(), sp(). Divides by alpha so the step is still 2841c alpha*ds(). 2842c 2843c This has little effect on most calculations but for (h2o)5 it 2844c reduces the number of iterations from 99 to 69. 2845c 2846 integer i 2847 logical gsopt_geom_cart_coords_get 2848c 2849 if (odebug) then 2850 write(6,*) ' Expected ds ' 2851 call doutput(ds, 1, nvar, 1, 1, nvar, 1, 1) 2852 endif 2853c 2854 call dcopy(nvar, sp, 1, ds, 1) ! Old coordinates into ds() 2855 if (.not. gsopt_geom_cart_coords_get(geom, sp)) 2856 $ call errquit('gsopt: geom?',0, GEOM_ERR) 2857c 2858 do i = 1, nvar 2859 ds(i) = sp(i) - ds(i) 2860 enddo 2861 if (odebug) then 2862 write(6,*) ' Actual ds ' 2863 call doutput(ds, 1, nvar, 1, 1, nvar, 1, 1) 2864 endif 2865c 2866 end 2867 logical function gsopt_geom_cart_coords_get(geom, coords) 2868 implicit none 2869#include "errquit.fh" 2870#include "geom.fh" 2871#include "nwc_const.fh" 2872#include "cgsopt.fh" 2873 integer geom 2874 double precision coords(*) 2875c 2876c If we are doing a periodic system and not using internals 2877c then we want the fractional coordinates. Otherwise cartesian. 2878c 2879c If we are including stress append the amatrix 2880c 2881 if (.not. geom_cart_coords_get(geom, coords)) 2882 $ call errquit('gsopt: geom cart?',0, GEOM_ERR) 2883c 2884 if (isystype.ne.0 .and. (.not. zcoord)) then 2885 if (.not. geom_cart_to_frac(geom, coords)) 2886 $ call errquit('gsopt: frac_to_cart?',0, GEOM_ERR) 2887 endif 2888c 2889 if (ostress) then 2890 if (.not. geom_amatrix_get(geom, coords(3*nat_real+1))) 2891 $ call errquit('gsopt: failed to get amatrix',0,0) 2892 endif 2893 if (ostress2) then 2894 if (.not. geom_lattice_get(geom, coords(3*nat_real+1))) 2895 $ call errquit('gsopt: failed to get lattice',0,0) 2896 endif 2897c 2898 gsopt_geom_cart_coords_get = .true. 2899c 2900 end 2901 logical function gsopt_geom_cart_coords_set(geom, coords) 2902 implicit none 2903#include "errquit.fh" 2904#include "geom.fh" 2905#include "nwc_const.fh" 2906#include "cgsopt.fh" 2907 integer geom 2908 double precision coords(*) 2909c 2910c If we are doing a periodic system and not using internals 2911c then we want the fractional coordinates. Otherwise cartesian. 2912c 2913 logical geom_amatrix_set 2914 external geom_amatrix_set 2915c 2916 if (ostress) then 2917 if (.not. geom_amatrix_set(geom, coords(3*nat_real+1))) 2918 $ call errquit('gsopt: failed to set amatrix',0,0) 2919 endif 2920 if (ostress2) then 2921 if (.not. geom_lattice_set(geom, coords(3*nat_real+1))) 2922 $ call errquit('gsopt: failed to set lattice',0,0) 2923 endif 2924 2925 if (isystype.ne.0 .and. (.not. zcoord)) then 2926 if (.not. geom_frac_to_cart(geom, coords)) 2927 $ call errquit('gsopt: frac_to_cart?',0,0) 2928 endif 2929 if (.not. geom_cart_coords_set(geom, coords)) 2930 $ call errquit('gsopt: geom cart?',0,0) 2931 if (isystype.ne.0 .and. (.not. zcoord)) then 2932 if (.not. geom_cart_to_frac(geom, coords)) 2933 $ call errquit('gsopt: frac_to_cart?',0,0) 2934 endif 2935c 2936 gsopt_geom_cart_coords_set = .true. 2937c 2938 end 2939 subroutine gsopt_get_grad(rtdb,geom) 2940 implicit none 2941#include "errquit.fh" 2942#include "mafdecls.fh" 2943#include "util.fh" 2944#include "nwc_const.fh" 2945#include "cgsopt.fh" 2946#include "rtdb.fh" 2947#include "geom.fh" 2948 integer rtdb, geom 2949 character*32 theory 2950c 2951c Get the gradient. 2952c 2953c If the optimization is supposed to be happening in fractional 2954c coordinates convert the gradients from cartesians. 2955c 2956c If we are including stress append the cell param gradients 2957c 2958 logical geom_grad_cart_to_frac 2959c 2960 if (.not. rtdb_get(rtdb, 'task:gradient', mt_dbl, ncart, 2961 $ gx)) call errquit('gsopt: could not get gradient',0,0) 2962 if (isystype .ne. 0) then 2963 if (.not. geom_grad_cart_to_frac(geom, gx)) 2964 $ call errquit('gsopt: frac_to_cart?',0,0) 2965 end if 2966 if (ostress) then 2967 if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory)) 2968 $ call errquit('gsopt: stress theory not specified',0,RTDB_ERR) 2969 if (theory.eq.'pspw') then 2970 if (.not. rtdb_get(rtdb, 'pspw:stress', mt_dbl, 9, 2971 $ gx(3*nat_real+1))) call errquit 2972 $ ('gsopt: could not get stress',0,0) 2973 else if (theory.eq.'band') then 2974 if (.not. rtdb_get(rtdb, 'band:stress', mt_dbl, 9, 2975 $ gx(3*nat_real+1))) call errquit 2976 $ ('gsopt: could not get stress',0,0) 2977 else if (theory.eq.'paw') then 2978 if (.not. rtdb_get(rtdb, 'paw:stress', mt_dbl, 9, 2979 $ gx(3*nat_real+1))) call errquit 2980 $ ('gsopt: could not get stress',0,0) 2981 else 2982 call errquit('gsopt: no stress in theory',0,RTDB_ERR) 2983 end if 2984 endif 2985 2986 if (ostress2) then 2987 if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory)) 2988 $ call errquit('gsopt: stress theory not specified',0,RTDB_ERR) 2989 if (theory.eq.'pspw') then 2990 if (.not. rtdb_get(rtdb, 'pspw:lstress', mt_dbl, 6, 2991 $ gx(3*nat_real+1))) call errquit 2992 $ ('gsopt: could not get stress',0,0) 2993 else if (theory.eq.'band') then 2994 if (.not. rtdb_get(rtdb, 'band:lstress', mt_dbl, 6, 2995 $ gx(3*nat_real+1))) call errquit 2996 $ ('gsopt: could not get stress',0,0) 2997 else if (theory.eq.'paw') then 2998 if (.not. rtdb_get(rtdb, 'paw:lstress', mt_dbl, 6, 2999 $ gx(3*nat_real+1))) call errquit 3000 $ ('gsopt: could not get stress',0,0) 3001 else 3002 call errquit('gsopt: no stress in theory',0,RTDB_ERR) 3003 end if 3004 endif 3005c 3006 end 3007cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3008cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3009 subroutine gsopt_pickstp(rtdb, geom, istep) 3010 implicit none 3011#include "errquit.fh" 3012#include "rtdb.fh" 3013#include "nwc_const.fh" 3014#include "cmepgs.fh" 3015#include "cgsopt.fh" 3016#include "mafdecls.fh" 3017#include "global.fh" 3018#include "util.fh" 3019#include "stdio.fh" 3020 integer rtdb 3021 integer geom 3022 integer istep 3023c 3024c this routine for minimization 3025c 3026c put into ds() a search direction in the optimization 3027c variables (internal or cartesian) based upon the 3028c current gradient, g(), and hessian. apply constraints. 3029c 3030c apply step restrictions by recommending an initial 3031c value for the line search parameter alpha. 3032c 3033c only alpha and ds() are modified. 3034c 3035 integer i, iat 3036cjmc 3037 integer iact, lowest 3038 double precision lambda, tau, term 3039 double precision tiniest 3040 double precision big 3041 parameter (tiniest = 1.0e-14) 3042 parameter (big = 1d6) 3043cjmc 3044* integer info 3045 integer l_hess, k_hess, l_work, k_work, lenwork 3046 double precision eigval(max_nvar) ! hessian eigenvalues 3047 double precision pg(max_nvar) ! p.g 3048 double precision pr(max_nvar) ! p.radius 3049 double precision gv(max_nvar) ! gradient along eigenvectors 3050 double precision dv(max_nvar) ! step along eigenvectors 3051 double precision gc(max_nvar) ! step along eigenvectors 3052 double precision coords(max_nvar) ! step along eigenvectors 3053 double precision dsmax ! max. value of current step (smax is prev.) 3054c 3055 double precision beta, s0g0, s0g1, s1g0, s1g1, numerator, 3056 $ denominator 3057 double precision bohr, deg ! for printing purposes 3058 double precision trustds ! restriction of step in opt. variable 3059 logical ophigh 3060 logical gsopt_geom_cart_coords_get 3061c 3062c get the hessian and gradient with appropriate projectors 3063c applied following peng, ayala, schlegel and frisch so that 3064c redundant internal modes are shifted to high eigenvalues. 3065c 3066 ophigh = util_print('high', print_high) 3067 if (.not. ma_push_get(mt_dbl, nvar**2, 'hess', 3068 $ l_hess, k_hess)) call errquit 3069 $ ('gsopt_pickstp: memory for hessian',nvar**2, ma_err) 3070 call gsopt_project_hess_grad(dbl_mb(k_hess), pg) 3071c 3072c diagonalize the hessian. should really do the generalized 3073c eigenvalue problem since the underlying basis is not independent 3074c (if we are using autoz). not yet being done. 3075c 3076c to cause degenerate eigenvalues to be resolved into symmetry 3077c adapted combinations use jacobi not dsyev and screen out junk 3078c 3079 lenwork = max(nvar**2,100) 3080 if (.not. ma_push_get(mt_dbl, lenwork, 'work', 3081 $ l_work, k_work)) call errquit 3082 $ ('gsopt_pickstp: memory for hessian', lenwork, ma_err) 3083 do i = 0, nvar**2-1 3084 if (abs(dbl_mb(k_hess+i)).lt.1d-8) dbl_mb(k_hess+i) = 0d0 3085 enddo 3086c 3087c have eigenvalues in eigval, eigenvectors in dbl_mb(k_hess). 3088c 3089 call util_jacobi(nvar, dbl_mb(k_hess), nvar, eigval) 3090 if (odebug .or. (util_print('hvecs',print_never) 3091 $ .and. ga_nodeid().eq.0)) then 3092 write(6,*) ' Eigenvalues of the hessian ' 3093 call doutput(eigval, 1, nvar, 1, 1, nvar, 1, 1) 3094 write(6,*) ' Eigenvectors of the hessian ' 3095 call output(dbl_mb(k_hess), 1, nvar, 1, nvar, nvar, nvar, 1) 3096 endif 3097c 3098c project the gradient onto the hessian eigenvectors 3099c 3100 call dgemv('t', nvar, nvar, 1d0, dbl_mb(k_hess), nvar, 3101 $ pg, 1, 0d0, gv, 1) 3102c 3103c *** calculate constrained saddle gradient *** 3104c 3105 call dfill(nvar, 0d0, pr, 1) 3106 call dfill(nvar, 0d0, radius, 1) 3107c 3108 if (.not. gsopt_geom_cart_coords_get(geom, coords)) 3109 $ call errquit('tropt: geom?',0, GEOM_ERR) 3110c 3111 call dcopy(nvar, coords, 1, radius, 1) 3112 call daxpy(nvar, -1.0d0, center, 1, radius, 1) 3113c 3114c 3115CCCCCCCCCCCCCCCCCCCCCC 3116CCCC MASS CCCC 3117CCCCCCCCCCCCCCCCCCCCCC 3118 if (mswg) call mwcoord(radius, nvar, .true.) 3119CCCCCCCCCCCCCCCCCCCCCC 3120CCCC MASS CCCC 3121CCCCCCCCCCCCCCCCCCCCCC 3122c 3123c 3124 call dgemv('t', nvar, nvar, 1d0, dbl_mb(k_hess), nvar, 3125 $ radius, 1, 0d0, pr, 1) 3126c 3127c *** calculate combined projected cartesian gradient *** 3128c 3129 do iact=1,nvar 3130 gc(iact) = eigval(iact)*pr(iact) - gv(iact) 3131 end do 3132! 3133! *** find lowest (countable) eigenvalue position *** 3134! *** assuming a non-shifted hessian matrix *** 3135! 3136 do iact=1,nvar 3137 if (abs(eigval(iact)).gt.1.0d-8) then 3138 lowest = iact 3139 goto 100 3140 end if 3141 end do 3142c 3143 100 continue 3144! 3145! *** test for the "hard case", i.e. f1 + b1*d1 = 0 *** 3146! 3147 if (abs(gc(lowest)).lt.1.e-8) then 3148! 3149! *** set lambda to lowest eigenvalue *** 3150! 3151 lambda = -eigval(lowest) 3152! 3153! *** calculate reduced constrained step *** 3154! 3155 call dfill(nvar, 0d0, dv, 1) 3156 do iact=1,nvar 3157 if (iact.ne.lowest) then 3158 dv(iact) = gv(iact) + lambda*pr(iact) 3159 dv(iact) = -dv(iact)/(eigval(iact) + lambda + tiniest) 3160 end if 3161 end do 3162! 3163! *** check step length and modify, if necessary *** 3164! 3165 term = ddot(nvar, dv, 1, dv, 1) 3166 if (term.le.trust**2) then 3167 tau = sqrt(trust**2 - term) 3168 dv(lowest) = tau 3169 if (oprint) write(6,*) "hard case encountered" 3170 slength = sqrt(ddot(nvar, dv, 1, dv, 1)) 3171 goto 200 3172 end if 3173 end if 3174! 3175 lambda = 0.0 3176! 3177! *** restrict step size *** 3178! 3179 call gsopt_lambda(gc, eigval, lambda) 3180 call dfill(nvar, 0d0, dv, 1) 3181 do iact=1,nvar 3182 dv(iact) = gv(iact) + lambda*pr(iact) 3183 dv(iact) = -dv(iact)/(eigval(iact) + lambda + tiniest) 3184 end do 3185 slength = sqrt(ddot(nvar, dv, 1, dv, 1)) 3186! 3187! *** print shift factor(s) *** 3188! 3189 200 continue 3190 3191 if (oprint) then 3192 if (ga_nodeid().eq.0) write (6,5030) lambda 3193 5030 format ('lambda for step: ',g10.3) 3194 if (ga_nodeid().eq.0) write (6,5040) slength 3195 5040 format ('step magnitude : ',g10.3) 3196 end if 3197 3198 if (odebug) then 3199 write(6,*) ' Step in spectral form ' 3200 call doutput(dv, 1, nvar, 1, 1, nvar, 1, 1) 3201 endif 3202c 3203c transform back to optimization space 3204c 3205 call dgemv('n', nvar, nvar, 1d0, dbl_mb(k_hess), nvar, 3206 $ dv, 1, 0d0, ds, 1) 3207 if (odebug) then 3208 write(6,*) ' Step in optimization variables' 3209 call doutput(ds, 1, nvar, 1, 1, nvar, 1, 1) 3210 endif 3211c 3212c enforce symmetry 3213c 3214 call gsopt_symmetrize_step(geom) 3215c 3216c enforce frozen atoms in cartesians 3217c 3218 if (ga_nodeid().eq.0.and.ophigh) 3219 $ write(6,*) 'Zeroing constrained gradient' 3220 if ((.not. zcoord) .and. (nactive .ne. nat_real)) then 3221 do iat = 1, nat 3222 if (.not. oactive(iat)) then 3223 do i = 1, 3 3224 ds((iat-1)*3+i) = 0.0 3225 end do 3226 end if 3227 end do 3228 end if 3229c 3230 if (.not. ma_chop_stack(l_hess)) call errquit 3231 $ ('gsopt_pickstp: ma?',0, MA_ERR) 3232c 3233c edo seems to have encountered a case where different processors 3234c generated different steps. to prevent this, broadcast the 3235c critical info to everyone. 3236c 3237 call ga_brdcst(1,ds,8*nvar,0) 3238 call ga_brdcst(2,alpha,8,0) 3239c 3240 if (util_print('searchdir',print_high) .and. 3241 $ ga_nodeid().eq.0) then 3242 write(6,*) 3243 write(6,*) ' the search direction' 3244 call output(ds,1,3,1,nat,3,nat,1) 3245 write(6,*) 3246 call util_flush(6) 3247 endif 3248c 3249c 3250CCCCCCCCCCCCCCCCCCCCCC 3251CCCC MASS CCCC 3252CCCCCCCCCCCCCCCCCCCCCC 3253 if (mswg) then 3254 call mwcoord(radius, nvar, .false.) 3255 call mwcoord( ds, nvar, .false.) 3256 end if 3257CCCCCCCCCCCCCCCCCCCCCC 3258CCCC MASS CCCC 3259CCCCCCCCCCCCCCCCCCCCCC 3260c 3261c 3262 end 3263cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3264cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3265 subroutine gsopt_raphson(de1,dr,eigval,lambda) 3266! do newton-raphson step with hessian eigenvalue shift. 3267 3268 implicit none 3269#include "nwc_const.fh" 3270#include "cgsopt.fh" 3271 3272 double precision eps 3273 parameter (eps = 1d-5) 3274 double precision tiniest 3275 parameter (tiniest = 1.0d-14) 3276 3277 double precision de1(nvar),dr(nvar),eigval(nvar) 3278 3279 integer iact 3280 double precision lambda,term 3281! 3282! *** calculate displacement vector *** 3283! 3284 call dfill(nvar, 0d0, dr, 1) 3285 3286 do iact=1,nvar 3287 dr(iact) = -de1(iact)/(eigval(iact) + lambda + tiniest) 3288 end do 3289 3290 end 3291cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3292cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3293 subroutine gsopt_lambda(f,eigval,lambda) 3294! calculate hessian eigenvalue shift factor from 3295! trust-region quadratic approximation. 3296 3297 implicit none 3298#include "errquit.fh" 3299#include "global.fh" 3300#include "mafdecls.fh" 3301#include "stdio.fh" 3302#include "util.fh" 3303#include "nwc_const.fh" 3304#include "cgsopt.fh" 3305#include "rtdb.fh" 3306#include "geom.fh" 3307 3308 integer maxiter 3309 parameter (maxiter = 100) 3310 3311 double precision big,eps 3312 parameter (big = 1d6, eps = 1d-10) 3313 3314 double precision f(nvar),eigval(nvar) 3315 3316 logical error 3317 integer iact,iter,loop 3318 double precision check,dfn,dx,fn,itlamb,lambda,lower,upper 3319 double precision tmp1,tmp2 3320! 3321! *** initialization *** 3322! 3323 error = .false. 3324 upper = big 3325c 3326c *** search interval definition *** 3327c 3328 do iact=1,nvar 3329 if (abs(eigval(iact)).gt.1.e-8) then 3330 lower = -eigval(iact) 3331 goto 100 3332 end if 3333 end do 3334c 3335 100 continue 3336c 3337 itlamb = lower + 0.05 3338! 3339! *** iteration loop *** 3340! 3341 iter = 0 3342 3343 do loop=1,maxiter 3344! 3345! *** calculate fn(lambda) function *** 3346! 3347 fn = 0.0 3348 do iact=1,nvar 3349 fn = fn + (f(iact)/(eigval(iact) + itlamb))**2 3350 end do 3351c 3352c *** derivative of fn(lambda) with respect to lambda *** 3353c 3354 dfn = 0.0 3355 do iact=1,nvar 3356 dfn = dfn + f(iact)**2/(eigval(iact) + itlamb)**3 3357 end do 3358! 3359! *** calculate next estimate of lambda *** 3360! 3361 dx = fn/dfn 3362 lambda = itlamb + (sqrt(fn)/trust - 1.0)*dx 3363! 3364! *** select new lambda value *** 3365! 3366 if ((lambda.le.upper).and.(lambda.gt.lower)) then 3367 3368 if (lambda.lt.itlamb) then 3369 upper = itlamb 3370 else 3371 lower = itlamb 3372 end if 3373 3374 itlamb = lambda 3375 3376 else 3377 3378 if (lambda.lt.itlamb) then 3379 upper = itlamb 3380 else 3381 lower = itlamb 3382 end if 3383 3384 if (upper.eq.big) then 3385 itlamb = lambda + 0.05 3386 else 3387 itlamb = 0.5*(upper + lower) 3388 end if 3389 3390 end if 3391 3392 iter = iter + 1 3393 check = abs(sqrt(fn) - trust) 3394 if (check.le.eps) go to 200 3395 3396 end do 3397! 3398! *** iteration failed *** 3399! 3400 error = .true. 3401 3402 if (oprint.and.ga_nodeid().eq.0) write (6,5000) 34035000 format ('lambda iteration did not converge') 3404 3405 200 continue 3406! 3407! *** failed lambda calculation *** 3408! 3409 if (error) then 3410 write (6,'(a,f16.8)') 'lambda', lambda 3411 write (6,'(a,f16.8)') 'trust radii', trust 3412 call errquit('gsopt_lambda: fatal error', 911, geom_err) 3413 end if 3414 3415 end 3416CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCccc 3417CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCccc 3418 logical function mepgs_freq(rtdb) 3419c 3420c **** A copy of task_freq **** 3421c 3422 implicit none 3423#include "errquit.fh" 3424#include "mafdecls.fh" 3425#include "rtdb.fh" 3426#include "geom.fh" 3427#include "stdio.fh" 3428#include "global.fh" 3429#include "inp.fh" 3430#include "util.fh" 3431c 3432 logical task_hessian 3433 external task_hessian 3434c 3435 integer rtdb 3436c 3437 integer nat, geom 3438 logical status 3439 character*(nw_max_path_len) filehess 3440c 3441 double precision cpu, wall 3442c 3443 logical ignore 3444 logical o_reuse 3445c 3446 call ecce_print_module_entry('task frequencies') 3447 cpu = util_cpusec() 3448 wall = util_wallsec() 3449c 3450 if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .false.)) 3451 $ call errquit('task_freq: failed to invalidate status',0, 3452 & RTDB_ERR) 3453 if (ga_nodeid().eq.0 .and. 3454 $ util_print('task_freq', print_low)) then 3455 write(LuOut,*) 3456 write(LuOut,*) 3457 call util_print_centered(6, 3458 $ 'NWChem Nuclear Hessian and Frequency Analysis', 3459 $ 40,.true.) 3460 write(LuOut,*) 3461 endif 3462* 3463 if (rtdb_get(rtdb,'vib:reuse',mt_log,1,ignore)) then 3464 o_reuse = ignore 3465 else 3466 o_reuse = .false. 3467 endif 3468* 3469 if (.not.(o_reuse)) then 3470 status = task_hessian(rtdb) 3471 if (.not.status) call errquit 3472 & ('task_freq: task_hessian failed',911, CALC_ERR) 3473 else 3474 if (ga_nodeid().eq.0) 3475 & call util_print_centered(LuOut, 3476 & 'reusing previously generated Hessian', 3477 & 40,.true.) 3478 status = .true. 3479 endif 3480* 3481 ignore = rtdb_parallel(.false.) 3482 if ((ga_nodeid()).eq.0) then 3483 if (o_reuse) then 3484 if (rtdb_cget(rtdb,'vib:reuse_hessian_file',1,filehess)) then 3485 write(LuOut,*)' re-using hessian in file ', 3486 & filehess(1:inp_strlen(filehess)) 3487 else 3488* in case of manual restart and renaming of hess file 3489 call util_file_name('hess', .false., .false.,filehess) 3490 endif 3491 else 3492 if (.not. rtdb_cget(rtdb, 'task:hessian file name', 1, 3493 $ filehess)) call errquit 3494 $ ('task_freq: failed reading hessian filename from rtdb',0, 3495 & RTDB_ERR) 3496 endif 3497c 3498c create/load reference geometry 3499c 3500 if (.not.geom_create(geom,'geometry')) call errquit 3501 $ ('task_freq:geom_create failed?',1, GEOM_ERR) 3502 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) 3503 $ call errquit 3504 $ ('task_freq:geom_rtdb_load failed?',2, GEOM_ERR) 3505 if (.not. geom_ncent(geom,nat)) call errquit 3506 $ ('task_freq:geom_ncent failed?',3, GEOM_ERR) 3507 if (.not. geom_destroy(geom)) call errquit 3508 $ ('task_freq:geom_destroy failed?',911, GEOM_ERR) 3509 call mepgs_vib(rtdb,filehess,.true., 3510 $ 0,.false.,0,.false.,nat) 3511 endif 3512 call ga_sync() 3513 ignore = rtdb_parallel(.true.) 3514c 3515 cpu = util_cpusec() - cpu 3516 wall = util_wallsec() - cpu 3517c 3518 if (.not. rtdb_put(rtdb, 'task:cputime', mt_dbl, 1, cpu)) 3519 $ call errquit('task_freq: failed storing cputime',0, RTDB_ERR) 3520 if (.not. rtdb_put(rtdb, 'task:walltime', mt_dbl, 1, wall)) 3521 $ call errquit('task_freq: failed storing walltime',0, 3522 & RTDB_ERR) 3523 if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .true.)) 3524 $ call errquit('task_freq: failed to set status',0, 3525 & RTDB_ERR) 3526c 3527c 3528 call ecce_print1('cpu time', mt_dbl, cpu, 1) 3529 call ecce_print1('wall time', mt_dbl, wall, 1) 3530 mepgs_freq = status 3531c 3532c 3533 end 3534CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc 3535CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcc 3536 subroutine mepgs_vib(rtdb,hess_file,in_file,hess_ma,in_ma, 3537 & hess_ga,in_ga,natomin) 3538c 3539c **** A copy of vib_vib routine ***** 3540c 3541 IMPLICIT NONE ! REAL*8 (A-H,O-Z) 3542#include "errquit.fh" 3543 LOGICAL PROJEC,ZEROPE,HESOUT,INTERN 3544* CHARACTER*7 INPFIL 3545 INTEGER NATOM, NAT3, NHESS, NHESST 3546 COMMON /cvib_HESS/ NATOM,NAT3,NHESS,NHESST ! hessian information 3547#include "stdio.fh" 3548#include "mafdecls.fh" 3549#include "rtdb.fh" 3550c:: passed 3551 integer rtdb ! [input] rtdb handle 3552 character*(*) hess_file ! [input] name of file storing lower triangular packed hessian 3553 integer hess_ma ! [input] MA handle to square hessian 3554 integer hess_ga ! [input] GA handle to square hessian 3555 logical in_file ! [input] hessian is in file get it there 3556 logical in_ma ! [input] hessian is in MA array 3557 logical in_ga ! [input] hessian is in GA array 3558 integer natomin ! [input] number of atoms 3559c 3560 logical status 3561c 3562 integer i_core, h_core, iii, ioldlabs, ivc, itot 3563 integer nels, npri, ihess, icoord, ihesst, ihesstcp, 3564 & ihessp, iegval, iegvec, iddpol, iddpolq, intense 3565 integer imass, iscr, i10, i20, i30, i40,i_w1,l_w1,i_w2,l_w2 3566 double precision dbl_tmp 3567 logical first_pass 3568 character*255 dipole_file 3569 logical dipole_file_exists 3570 logical animation_on 3571c 3572 first_pass = .true. 3573*... check input logic 3574 status = (in_file.and.(in_ma.or.in_ga)) 3575 status = status.or.(in_ma.and.(in_file.or.in_ga)) 3576 status = status.or.(in_ga.and.(in_file.or.in_ma)) 3577 if (status) then 3578 write(luout,*)' ERROR: more than one source for hessian ' 3579 write(luout,*)' in_file :',in_file 3580 write(luout,*)' in_ma :',in_ma 3581 write(luout,*)' in_ga :',in_ga 3582 call errquit(' vib_vib: error ',911, UNKNOWN_ERR) 3583 endif 3584 if (in_ga) 3585 & call errquit 3586 & ('vib_vib: ga access to hessian not implemented yet',911, 3587 & CAPMIS_ERR) 3588C 3589C Zero core 3590C 3591 call vib_setup ! subroutine to set up some constants 3592 NATOM = natomin ! number of atoms in species. 3593 IF (NATOM.LE.1) THEN ! check for incorrect number of atoms 3594 WRITE(6,*)' You want to calculate the vibrational ', 3595 + 'frequencies for ',NATOM,' atoms?' 3596 WRITE(6,*)' Unfortunately this is not possible ' 3597 CALL errquit('vib_vib: bomb',911, INPUT_ERR) 3598 ENDIF 3599 NAT3 = NATOM*3 ! 3-N (as in degrees of freedom) 3600 NHESS = NAT3*NAT3 ! dimension of hessian 3601 NHESST = NAT3*(NAT3+1)/2 ! dimension of lower triangular hessian 3602 NELS = 7*MAX(3*NATOM-6,1) 3603 NPRI = 0 3604C 3605C Calculate pointers 3606C 3607 IHESS = 1 ! square hessian 3608 IHESST = IHESS + NHESS ! lower-tri Hessian 3609 ihesstcp = IHESST + NHESST ! copy of lower-tri hessian 3610 ICOORD = IHESSTcp + NHESST ! geometrical coordinates 3611 IMASS = ICOORD + NAT3 ! mass of each atom 3612 IEGVAL = IMASS + NATOM ! eigenvalues from Hessian matrix 3613 IEGVEC = IEGVAL + NAT3 ! eigenvectors from Hessian matrix 3614 ISCR = IEGVEC + NHESS ! dynamic bottom of core array 3615 IHESSP = ISCR + 8*NAT3 ! addition of scratch space needed 3616C--------The following are pointers for GAMESS internal coordinate subroutines. 3617C 3618 I10 = IHESSP + NAT3*NAT3 ! space for zmat 3619 I20 = I10 + NELS 3620 I30 = I20 + NAT3*NAT3 ! Space to represent internal coord. Hessian 3621 I40 = I30 + NAT3*NAT3 ! 3622 Iddpol = I40 + 8*NAT3 ! derivative dipole in cartesians 3623 Iddpolq = Iddpol + 3*NAT3 ! derivative dipole in normal modes 3624 Intense = Iddpolq+ 3*NAT3 ! intensities 3625 ITOT = Intense + 3*natom*4 3626 itot = itot + 2*natom+1 + 6*nat3 ! extra for call to rdinp 3627c 3628 if (.not.ma_push_get 3629 & (MT_DBL,itot,' core for vib ',h_core, i_core)) 3630 & call errquit('vib_vib: ma_push_get failed ',911, MA_ERR) 3631C 3632C Reset pointers for MA array 3633C 3634 IHESS = i_core ! square hessian 3635 IHESST = IHESS + NHESS ! lower-tri Hessian 3636 ihesstcp = IHESST + NHESST ! copy lower-tri Hessian 3637 ICOORD = IHESSTcp + NHESST ! geometrical coordinates 3638 IMASS = ICOORD + NAT3 ! mass of each atom 3639 IEGVAL = IMASS + NATOM ! eigenvalues from Hessian matrix 3640 IEGVEC = IEGVAL + NAT3 ! eigenvectors from Hessian matrix 3641 ISCR = IEGVEC + NHESS ! dynamic bottom of core array 3642 IHESSP = ISCR + 8*NAT3 ! addition of scratch space needed 3643C--------The following are pointers for GAMESS internal coordinate subroutines. 3644C 3645 I10 = IHESSP + NAT3*NAT3 ! space for zmat 3646 I20 = I10 + NELS 3647 I30 = I20 + NAT3*NAT3 ! Space to represent internal coord. Hessian 3648 I40 = I30 + NAT3*NAT3 ! 3649 ioldlabs = I40 + 8*NAT3 ! 3650 ivc = ioldlabs + 2*natom + 1 3651 iddpol = ivc + 6*nat3 3652 iddpolq = iddpol + 3*nat3 3653 intense = iddpolq + 3*nat3 3654 Itot = intense + 3*natom*4 3655c 3656c read/load hessian and form triangle/square as needed 3657c 3658 if (in_ma) then 3659 ihess = hess_ma ! simply reset ptr to dbl_mb 3660* form triangle 3661 call vib_dtrngl(dbl_mb(ihess),dbl_mb(ihesst),nat3,nat3) 3662 endif 3663 if (in_file) then 3664 open(unit=69,file=hess_file,form='formatted',status='old', 3665 & err=99900,access='sequential') 3666 do iii = 0,(nhesst-1) 3667 read(69,*,err=99901,end=99902)dbl_tmp 3668 dbl_mb(ihesst+iii) = dbl_tmp 3669 enddo 3670 close(unit=69,status='keep') 3671 call vib_dsquar(dbl_mb(ihesst),dbl_mb(ihess),nat3,nat3) 3672 endif 3673 call util_file_name('fd_ddipole',.false., .false.,dipole_file) 3674 dipole_file_exists = .false. 3675 inquire(file=dipole_file,exist=dipole_file_exists) 3676 if (dipole_file_exists) then 3677 open(unit=70,file=dipole_file,form='formatted',status='old', 3678 & err=89900,access='sequential') 3679 do iii = 0,((3*nat3)-1) 3680 read(70,*,err=89901,end=89902) dbl_tmp 3681 dbl_mb(iddpol+iii) = dbl_tmp 3682 enddo 3683 close(unit=70,status='keep') 3684 endif 368500001 continue 3686 write(luout,*) 3687 write(luout,*) 3688 if (.not. first_pass) then 3689 WRITE(luout,*) 3690 & ' Vibrational analysis via the FX method ' 3691 write(luout,*) 3692 & ' --- with translations and rotations projected out ---' 3693 write(luout,*) 3694 & ' --- via the Eckart algorithm ---' 3695 endif 3696 if (first_pass) then 3697c 3698c save a copy of hesst 3699c 3700 call dcopy(nhesst,dbl_mb(ihesst),1,dbl_mb(ihesstcp),1) 3701 else 3702c 3703c restore copy of hesst and hess 3704c 3705 call dcopy(nhesst,dbl_mb(ihesstcp),1,dbl_mb(ihesst),1) 3706 call vib_dsquar(dbl_mb(ihesst),dbl_mb(ihess),nat3,nat3) 3707 endif 3708C 3709C Read in user input and tape10 arrays. NO INPUT REQUIRED NOW 3710C Note: ! scratch pointer for atom charges (real NAT words) 3711C and atom lables (real 2*nat words) 3712 call vib_rdinp( 3713 & dbl_mb(ihess),dbl_mb(ihesst),dbl_mb(icoord), 3714 & dbl_mb(imass),dbl_mb(iscr), dbl_mb(ioldlabs), 3715 & dbl_mb(i10),nels,projec,zerope,hesout,intern, 3716 & rtdb,first_pass) 3717 if (projec) then 3718 if (.not.ma_push_get 3719 & (MT_DBL,nat3*nat3,' w1 ',l_w1, i_w1)) 3720 & call errquit('vib_vib: ma_push_get failed ',911, MA_ERR) 3721 if (.not.ma_push_get 3722 & (MT_DBL,nat3*nat3,' w2 ',l_w2, i_w2)) 3723 & call errquit('vib_vib: ma_push_get failed ',911, MA_ERR) 3724 call vib_eckart( dbl_mb(ihess), dbl_mb(ihessp), dbl_mb(ihesst), 3725 & dbl_mb(icoord), dbl_mb(ivc), dbl_mb(i_w1),dbl_mb(i_w2)) 3726 if(.not.ma_chop_stack(l_w1)) call errquit( 3727 ' ' vib_vib: machopstack failed',1, MA_ERR) 3728 end if 3729* rak dfill 3730 CALL Dfill(NAT3,0.0d00,DBL_MB(ISCR),1) ! zero scratch used 3731C 3732 CALL mepgs_vib_hmass(DBL_MB(IHESST),DBL_MB(IMASS)) ! mass weight and scale hessian 3733C 3734C Diagonalize mass-weighted, scaled hessian matrix 3735C Note: ! scratch pointer for givens (real 5*NAT3 words) 3736c use hessp as scratch now calling rsg 3737C 3738 CALL vib_CALLG(DBL_MB(IHESSt),nhesst,DBL_MB(IHESSP), 3739 & dbl_mb(iscr),dbl_mb(iscr+nat3),DBL_MB(IEGVAL), 3740 & DBL_MB(IEGVEC), NAT3,NAT3) 3741C 3742C 3743C 3744CJMC revisar 3745 if (first_pass) CALL vib_NMASS(DBL_MB(IEGVEC),DBL_MB(IMASS)) ! "unmass" weight the normal modes. 3746CJMC revisar 3747C 3748C *** Note: DBL_MB(IHESST) now destroyed if needed reinitialize from DBL_MB(IHESS) 3749C 3750* rak dfill 3751 call dfill(5*nat3,0.0d00,dbl_mb(iscr),1) ! zero scratch used 3752C 3753CJMC CALL vib_WRTFREQ(rtdb,DBL_MB(IEGVAL),NAT3,ZEROPE,NPRI) ! Write out the zero-point energy 3754C 3755 CALL vib_CLEAN(DBL_MB(IEGVEC),NAT3*NAT3,1.0D-27) ! CLEAN eigenvectors 3756 if (.not. first_pass .and .dipole_file_exists) then 3757 call vib_intense(rtdb,dbl_mb(iegvec),dbl_mb(iegval),natom, 3758 & dbl_mb(iddpol),dbl_mb(iddpolq),dbl_mb(intense), 3759 & first_pass) 3760 endif 3761c 3762 if (.not.first_pass) then 3763* if any negative eigenvalues print out steps in their direction 3764 call mepgs_vib_istep( 3765 & rtdb,nat3,natom, 3766 & dbl_mb(iegvec),dbl_mb(iegval), 3767 & dbl_mb(icoord),dbl_mb(iscr), 3768 & dbl_mb(iscr+nat3),dbl_mb(iscr+(2*nat3))) 3769 3770 endif 3771c 3772 if (first_pass) then 3773 first_pass = .false. 3774 goto 00001 3775 endif 3776C 3777 if (.not.ma_pop_stack(h_core)) call errquit 3778 & ('vib_rdinp ma_pop failed',911, MA_ERR) 3779 return 378089900 continue 3781 write(luout,*)'dipole_file => ',dipole_file 3782 call errquit('vib_vib: error opening file: "dipole_file"',811, 3783 & DISK_ERR) 378489901 continue 3785 write(luout,*)'dipole_file => ',dipole_file 3786 call errquit('vib_vib: error reading file: "dipole_file"',811, 3787 & DISK_ERR) 378889902 continue 3789 write(luout,*)'dipole_file => ',dipole_file 3790 call errquit 3791 & ('vib_vib: unexpected EOF when reading file: "dipole_file"',811, 3792 & DISK_ERR) 379399900 continue 3794 write(luout,*)'hess_file => ',hess_file 3795 call errquit('vib_vib: error opening file: "hess_file"',911, 3796 & DISK_ERR) 379799901 continue 3798 write(luout,*)'hess_file => ',hess_file 3799 call errquit('vib_vib: error reading file: "hess_file"',911, 3800 & DISK_ERR) 380199902 continue 3802 write(luout,*)'hess_file => ',hess_file 3803 call errquit 3804 & ('vib_vib: unexpected EOF when reading file: "hess_file"',911, 3805 & DISK_ERR) 3806 END 3807CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc 3808CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc 3809 subroutine mepgs_vib_istep(rtdb,nat3,natom, 3810 & eigenvecs,eigenvals,coords,steps,rawstep,master) 3811c 3812c **** Copy of vib_istep 3813c 3814 implicit none 3815#include "errquit.fh" 3816#include "stdio.fh" 3817#include "geom.fh" 3818#include "nwc_const.fh" 3819#include "cmepgs.fh" 3820 double precision ddot 3821 external ddot 3822* 3823 integer rtdb ! [input] rtdb handle 3824 integer natom ! [input] number of atoms 3825 integer nat3 ! [input] 3*number of atoms 3826 double precision eigenvecs(nat3,nat3) ! [input](xyz&atom,mode) 3827 double precision eigenvals(nat3) ! [input] (mode) 3828 double precision master(3,natom) ! [scratch] original coordintates 3829 double precision coords(3,natom) ! [scratch] coords after step 3830 double precision rawstep(3,natom) ! [scratch] step generated by vector 3831 double precision steps(3,natom) ! [scratch] step generated by vector 3832c 3833 integer imode,ivec,iatom,ixyz 3834 integer geom 3835 double precision scale 3836 double precision xyz(3),charge 3837 double precision length_of_step 3838 character*16 tag 3839 character*10 units 3840 intrinsic sqrt 3841c 3842CJMC 3843 double precision delta 3844 double precision atmass(natom) 3845CJMC 3846 double precision thresh 3847 parameter (thresh=1.0d-2) 3848c::-statement function 3849 logical is_it_close_to 3850 double precision value,test 3851 intrinsic abs 3852*--- is value close to test? 3853 is_it_close_to(value,test) = (abs(value-test).lt.thresh) 3854c 3855 if (.not.geom_create(geom,'geometry')) call errquit 3856 & ('vib_istep: geom create failed',911, GEOM_ERR) 3857 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit 3858 & ('vib_istep: geom_rtdb_load failed',911, RTDB_ERR) 3859 if (.not.geom_cart_coords_get(geom,master)) call errquit 3860 & ('vib_istep: geom_get_cart_coords failed',911, GEOM_ERR) 3861 if (.not.geom_get_user_scale(geom,scale)) call errquit 3862 & ('vib_istep: geom_get_user_scale failed',911, GEOM_ERR) 3863 if (.not.geom_get_user_units(geom,units)) call errquit 3864 & ('vib_istep: geom_get_user_units failed',911, GEOM_ERR) 3865c 3866 imode = 1 3867 ivec = 0 3868 write(luout,10000)imode,eigenvals(imode) 3869 call dfill(nat3,0.0d00,rawstep,1) 3870 do iatom = 1, natom 3871 do ixyz = 1,3 3872 ivec = ivec+1 3873 rawstep(ixyz,iatom) = eigenvecs(ivec,imode) 3874 enddo 3875 enddo 3876CJMC 3877 delta = sqrt(2.0d0*evib/abs(eigenvals(imode))) 3878CJMC 3879 call dcopy(nat3,rawstep,1,steps,1) 3880 call dscal(nat3, delta, steps, 1) 3881 length_of_step = sqrt(ddot(nat3,steps,1,steps,1)) 3882 write(luout,10001)length_of_step,units 3883 do iatom=1,natom 3884 if(.not.geom_cent_get(geom,iatom,tag,xyz,charge)) 3885 & call errquit('vib_istep: geom_cent_get failed',911, GEOM_ERR) 3886 write(luout,10002)iatom,tag,charge, 3887 & (steps(ixyz,iatom),ixyz=1,3) 3888 enddo 3889 write(luout,10003) 3890c 3891c **** Mass weight coordinates **** 3892c 3893 if (.not. geom_masses_get(geom,natom,atmass)) 3894 & call errquit('vib_step: geom_masses_get failed',911, GEOM_ERR) 3895 do iatom=1, natom 3896 do ixyz=1,3 3897 master(ixyz,iatom) = master(ixyz,iatom)*sqrt(atmass(iatom)) 3898 end do 3899 end do 3900c 3901c **** Store "forward" direction **** 3902c 3903 call dcopy(nat3,master,1,coords,1) 3904 call daxpy(nat3,1.0d00,steps,1,coords,1) 3905c 3906c **** Unmass weight coordinates **** 3907c 3908 do iatom=1, natom 3909 do ixyz=1,3 3910 coords(ixyz,iatom) = coords(ixyz,iatom)/sqrt(atmass(iatom)) 3911 end do 3912 end do 3913 if (.not.geom_cart_coords_set(geom,coords)) call errquit 3914 & ('vib_istep: geom_cart_coords_set failed',911, GEOM_ERR) 3915 if (.not. geom_rtdb_store(rtdb, geom, 'ircforward')) 3916 $ call errquit('mepgs_vib_istep: grs?',geom, RTDB_ERR) 3917c 3918c **** Store "backward" direction **** 3919c 3920 call dcopy(nat3,master,1,coords,1) 3921 call daxpy(nat3,-1.0d00,steps,1,coords,1) 3922c 3923c **** Unmass weight coordinates **** 3924c 3925 do iatom=1, natom 3926 do ixyz=1,3 3927 coords(ixyz,iatom) = coords(ixyz,iatom)/sqrt(atmass(iatom)) 3928 end do 3929 end do 3930 if (.not.geom_cart_coords_set(geom,coords)) call errquit 3931 & ('vib_istep: geom_cart_coords_set failed',911, GEOM_ERR) 3932 if (.not. geom_rtdb_store(rtdb, geom, 'ircbackward')) 3933 $ call errquit('mepgs_vib_istep: grs?',geom, RTDB_ERR) 3934c 3935c **** Deallocate geom **** 3936c 3937 if (.not.geom_destroy(geom)) call errquit 3938 & ('vib_istep: geom_destroy failed',911, GEOM_ERR) 3939c 3940c 394110000 format(/,/,/,1x,78('='),/,6x,'Negative Nuclear Hessian Mode', 3942 & i5,2x,'Eigenvalue = ',f9.2,' a.u.',/,1x,78('-')) 394310001 format(2x,' Raw step length:',f7.3,1x,a10,';', 3944 & 2x, 'The Raw step for this mode is:') 394510002 format(' ',i4,' ',a16,' ',f10.4,3f15.8) 394610003 format(78('-')) 3947 end 3948CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3949CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3950 SUBROUTINE mepgs_vib_HMASS(HESST,ATMASS) 3951* $Id$ 3952C 3953C This routine mass weights and scales the Hessian matrix 3954C The scaling is done to avoid numerical problems in the 3955C diagonalization routine 3956C 3957 IMPLICIT NONE ! REAL*8 (A-H,O-Z) 3958#include "cmepgs.fh" 3959 INTEGER NAT, NAT3, NHESS, NHESST 3960 COMMON /CVIB_HESS/ NAT,NAT3,NHESS,NHESST ! HESSIAN INFORMATION 3961c 3962 double precision HESST(NHESST) ! lower triangular Hessian 3963 double precision ATMASS(NAT) ! mass of the atoms 3964CCCCCCCCCCCCCCCCCCCCCC 3965CCCC MASS CCCC 3966CCCCCCCCCCCCCCCCCCCCCC 3967 double precision mwhess(nat3*nat3) 3968CCCCCCCCCCCCCCCCCCCCCC 3969CCCC MASS CCCC 3970CCCCCCCCCCCCCCCCCCCCCC 3971c 3972 double precision fact, scale 3973 integer ii, jj, jjend, iatii, iatjj, idum 3974 double precision mass_ii, mass_jj 3975C 3976C set up function for locating i,j elements packed canonically as ij 3977C 3978 integer i, j, isym2, iatom 3979 ISYM2(I,J)=MAX(I,J)*((MAX(I,J))-1)/2 + MIN(I,J) 3980 IATOM(I) = (I+2)/3 ! function for coordinate I is on atom IATOM 3981C 3982 DO 00100 II = 1,NAT3 ! loop over coordinates 3983 JJEND = II 3984 IATII = IATOM(II) ! coordinate II is for atom IATII 3985 DO 00200 JJ = 1,JJEND ! loop over coordinates 3986 IDUM = ISYM2(II,JJ) ! get canonical index 3987 IATJJ = IATOM(JJ) ! coordinate JJ is for atom IATJJ 3988 mass_ii = atmass(iatii) 3989 mass_jj = atmass(iatjj) 3990 if (abs(mass_ii).lt.1.0d-01) mass_ii = 1.0d05 3991 if (abs(mass_jj).lt.1.0d-01) mass_jj = 1.0d05 3992* FACT = SQRT(ATMASS(IATII))*SQRT(ATMASS(IATJJ)) ! mass weight 3993 FACT = SQRT(mass_ii)*SQRT(mass_jj) ! mass weight 3994 HESST(IDUM) = HESST(IDUM)/FACT ! weight Hessian 399500200 CONTINUE 3996* idum = isym2(ii,ii) ! get canonical index for diagonal element 3997* if (abs(hesst(idum)).lt.1.0d-10) hesst(idum) = 1.0d04 399800100 CONTINUE 3999CJMC 4000c *** Avoid scaling --> Gives better agreement *** 4001c *** with quadratic model for moving away from TS *** 4002CJMC 4003 SCALE = 1.D0 ! Hessian scaling factor 4004*dscal 4005 call dscal(nhesst,scale,hesst,1) ! Scale Hessian for diagonaization 4006 4007CCCCCCCCCCCCCCCCCCCCCC 4008CCCC MASS CCCC 4009CCCCCCCCCCCCCCCCCCCCCC 4010 call vib_dsquar(hesst, mwhess, nat3, nat3) 4011 call geom_hnd_put_data('irc.hess', mwhess, nat3**2) 4012CCCCCCCCCCCCCCCCCCCCCC 4013CCCC MASS CCCC 4014CCCCCCCCCCCCCCCCCCCCCC 4015 4016 END 4017CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4018CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4019 subroutine mwcoord(vector, ndeg, put_mass) 4020 implicit none 4021#include "errquit.fh" 4022#include "util.fh" 4023#include "nwc_const.fh" 4024#include "cmepgs.fh" 4025#include "cgsopt.fh" 4026 integer ndeg 4027 double precision vector(ndeg) 4028 logical put_mass 4029c 4030 integer ipos, iatom, ixyz 4031c 4032 if (put_mass) then 4033 ipos = 0 4034 do iatom=1, nat 4035 do ixyz=1,3 4036 ipos = ipos + 1 4037 vector(ipos) = vector(ipos)*sqrt(atmass(iatom)) 4038 end do 4039 end do 4040 else if (.not. put_mass) then 4041 ipos = 0 4042 do iatom=1, nat 4043 do ixyz=1,3 4044 ipos = ipos + 1 4045 vector(ipos) = vector(ipos)/sqrt(atmass(iatom)) 4046 end do 4047 end do 4048 end if 4049c 4050 end 4051 subroutine mwgrad(vector, ndeg, put_mass) 4052 implicit none 4053#include "errquit.fh" 4054#include "util.fh" 4055#include "nwc_const.fh" 4056#include "cmepgs.fh" 4057#include "cgsopt.fh" 4058 integer ndeg 4059 double precision vector(ndeg) 4060 logical put_mass 4061c 4062 integer ipos, iatom, ixyz 4063c 4064 if (put_mass) then 4065 ipos = 0 4066 do iatom=1, nat 4067 do ixyz=1,3 4068 ipos = ipos + 1 4069 vector(ipos) = vector(ipos)/sqrt(atmass(iatom)) 4070 end do 4071 end do 4072 else if (.not. put_mass) then 4073 ipos = 0 4074 do iatom=1, nat 4075 do ixyz=1,3 4076 ipos = ipos + 1 4077 vector(ipos) = vector(ipos)*sqrt(atmass(iatom)) 4078 end do 4079 end do 4080 end if 4081c 4082 end 4083 4084