1c $Id$ 2 subroutine vscfm(NGRID,DMDR,NCOUP,VCFCT,IEXC, 3 * rtdb,nat,vec,eigen) 4 implicit double precision(a-h,o-z) 5#include "errquit.fh" 6c 7#include "nwc_const.fh" 8#include "geomP.fh" 9#include "geom.fh" 10#include "global.fh" 11#include "mafdecls.fh" 12#include "stdio.fh" 13#include "rtdb.fh" 14c 15 character*255 dipole_file 16 logical dmdr, dipole_file_exists 17 logical linear, restart 18 integer geom, rtdb, nat, geom_ref, nnm, geom_new 19 dimension vec(3*nat,3*nat), eigen(3*nat) 20 integer NSTART(7) 21c 22c ----- allocate memory and carry out vscf ----- 23c 24 nc1 = 3*nat 25 nc2 = (nc1*nc1+nc1)/2 26 nc3 = nc1*nc1 27c 28c nnm2 is number of pairs of modes (not triangular storage) 29c 30 if (.not.rtdb_get(rtdb,'vib:linear',mt_log,1,linear)) 31 $ linear = .false. 32 nnm = 3*nat - 6 33 if(linear) nnm = 3*nat-5 34 nnm2 = (nnm*nnm-nnm)/2 35 nnm3 = nnm*nnm 36c 37c when we have only one frequency we can only do ncoup=1 38c 39 IF (NNM.eq.1) NCOUP=1 40c 41 ngrid2 = (ngrid*ngrid+ngrid)/2 42 ngrid3 = ngrid*ngrid 43c 44c number of vibr. states for which wavefunctions are stored 45c 46 NST=NGRID/2 47 IF(NST.LT.4) NST=4 48c 49c maximum excitation level used in vibrational MP2 50c 51 NMAX=NST-1 52c 53c number of virtual states in MP2 54c 55 NVIRST=NMAX*NMAX*NNM2 + NMAX*NNM + 1 56C 57C number of triples of normal modes 58C no. of virtual states in MP2 includes triple excitations 59c 60 if (ncoup.gt.2) then 61 NTR=NNM*(NNM-1)*(NNM-2)/6 62 NVIRST=NVIRST+NMAX*NMAX*NMAX*NTR 63 end if 64c 65c memory for vscf 66c the first two of these hold harmonic normal modes and frequencies 67c 68 if (.not. 69 & ma_push_get(mt_dbl,NNM*NNM*NNM*NGRID*NGRID*NGRID, 70 & 'XTRIPV',l_XTRIPV,k_XTRIPV)) 71 & call errquit('vscfm: could not allocate l_XTRIPV',555, MA_ERR) 72 if (.not. 73 & ma_push_get(mt_dbl,nnm,'vscf freqs',l_fvscf,k_fvscf)) 74 & call errquit('vscfm: could not allocate l_fvscf',555, MA_ERR) 75 if (.not. 76 & ma_push_get(mt_dbl,nnm,'mp2 freqs',l_fmp2,k_fmp2)) 77 & call errquit('vscfm: could not allocate l_fmp2',555, MA_ERR) 78 if (.not. 79 & ma_push_get(mt_dbl,nnm,'delta Qs',l_dq,k_dq)) 80 & call errquit('vscfm: could not allocate l_dq',555, MA_ERR) 81 if (.not. 82 & ma_push_get(mt_dbl,nnm*ngrid,'grid Qs',l_rq,k_rq)) 83 & call errquit('vscfm: could not allocate l_rq',555, MA_ERR) 84 if (.not. 85 & ma_push_get(mt_dbl,nnm*(nnm+1)*ngrid,'awave',l_awave,k_awave)) 86 & call errquit('vscfm: could not allocate l_awave',555, MA_ERR) 87c 88c Get dipole information from previous frequency calculation 89c 90 call util_file_name('fd_ddipole',.false., .false.,dipole_file) 91 dipole_file_exists = .false. 92 inquire(file=dipole_file,exist=dipole_file_exists) 93 if (DMDR.and.iexc.eq.1.and.dipole_file_exists) then 94 ndipmo=0 95 if (.not. 96 & ma_push_get(mt_dbl,3*nc1,'dipole deriv',l_ddm,k_ddm)) 97 & call errquit('vscfm: could not allocate l_ddm',555, MA_ERR) 98 open(unit=70,file=dipole_file,form='formatted',status='old', 99 & err=89900,access='sequential') 100 do iii = 0,((3*nc1)-1) 101 read(70,*,err=89901,end=89902) dbl_tmp 102 dbl_mb(k_ddm) = dbl_tmp 103 enddo 104 close(unit=70,status='keep') 105 else 106 ndipmo=1 107 if (.not. 108 & ma_push_get(mt_dbl,nnm*ngrid,'diag dipole x',l_dmx,k_dmx)) 109 & call errquit('vscfm: could not allocate l_dmx',555, MA_ERR) 110 if (.not. 111 & ma_push_get(mt_dbl,nnm*ngrid,'diag dipole y',l_dmy,k_dmy)) 112 & call errquit('vscfm: could not allocate l_dmy',555, MA_ERR) 113 if (.not. 114 & ma_push_get(mt_dbl,nnm*ngrid,'diag dipole z',l_dmz,k_dmz)) 115 & call errquit('vscfm: could not allocate l_dmz',555, MA_ERR) 116c 117c Calculation of dipole should be at center of mass 118c 119 if (.not.rtdb_put(rtdb,'prop:center',mt_int,1,2)) 120 & call errquit('vscfm: rtdb_put of dipole:center failed',555, 121 & RTDB_ERR) 122 endif 123c 124c Allocate additional memory 125c 126 if (.not. 127 & ma_push_get(mt_dbl,nnm,'scaled freqs',l_sfreq,k_sfreq)) 128 & call errquit('vscfm: could not allocate l_sfreq',555, MA_ERR) 129 if (.not. 130 & ma_push_get(mt_dbl,nnm*ngrid,'diagonal V',l_diagv,k_diagv)) 131 & call errquit('vscfm: could not allocate l_diagv',555, MA_ERR) 132 if (.not. 133 & ma_push_get(mt_dbl,nnm3*ngrid3,'coupled V',l_coupv,k_coupv)) 134 & call errquit('vscfm: could not allocate l_coupv',555, MA_ERR) 135 if (.not. 136 & ma_push_get(mt_dbl,nnm3*ngrid3,'coup dipole x',l_dm2x,k_dm2x)) 137 & call errquit('vscfm: could not allocate l_dm2x',555, MA_ERR) 138 if (.not. 139 & ma_push_get(mt_dbl,nnm3*ngrid3,'coup dipole y',l_dm2y,k_dm2y)) 140 & call errquit('vscfm: could not allocate l_dm2y',555, MA_ERR) 141 if (.not. 142 & ma_push_get(mt_dbl,nnm3*ngrid3,'coup dipole z',l_dm2z,k_dm2z)) 143 & call errquit('vscfm: could not allocate l_dm2z',555, MA_ERR) 144 if (.not. 145 & ma_push_get(mt_dbl,nc1,'delta Xs',l_dx,k_dx)) 146 & call errquit('vscfm: could not allocate l_dx',555, MA_ERR) 147c 148c Check if we can restart: 149c - check where where we are on diag when ncoup = 1, 2, 3 150c - check where where we are on coupl when ncoup = 1, 2 151c - check where where we are on tripl when ncoup = 3 152c - get reference geometry from rtdb 153c - if any of the diag, coupl, or tripl is incomplete, deterimine 154c where we were and start from that point, minimizing the overhead 155c 156c nstart will be used to house restart data: 157c nstart(1) = level that was completed, i.e. none=0, diag=1, coupl=2, tripl=3 158c nstart(2-7) = restart point in level, needed are mode and grid point 159c i.e., two values for diag, four for coupl and six for tripl 160c 161c if restart, get reference geometry 162c if not a restart, get the geometry and make a copy to save 163c while we do potential scans 164c 165 call vscf_restart(rtdb,restart,ncoup,nstart,ldiag,lcoup,ltrip, 166 & nnm,ngrid,dbl_mb(k_rq),dbl_mb(k_dq), 167 & dbl_mb(k_diagv),dbl_mb(k_dmx),dbl_mb(k_dmy), 168 & dbl_mb(k_dmz),dbl_mb(k_dm2x),dbl_mb(k_dm2y), 169 & dbl_mb(k_dm2z),dbl_mb(k_coupv), 170 & dbl_mb(k_XTRIPV),ndipmo) 171 if (restart) then 172 if (.not.geom_create(geom_ref,'vscf_reference')) call errquit 173 * ('vscfm: vscf_reference geom_create failed',551, GEOM_ERR) 174 if (.not.geom_rtdb_load(rtdb,geom_ref,'vscf_reference')) 175 * call errquit('vscfm: geom_rtdb_load failed',551, RTDB_ERR) 176 if (.not.geom_rtdb_store(rtdb,geom_ref,'geometry')) 177 * call errquit('vscfm: geom_rtdb_store failed',551, RTDB_ERR) 178 else 179 if (.not.geom_create(geom_ref,'geometry')) call errquit 180 * ('vscfm: geometry geom_create failed',552, GEOM_ERR) 181 if (.not.geom_rtdb_load(rtdb,geom_ref,'geometry')) 182 * call errquit('vscfm: geometry geom_rtdb_load failed',552, 183 & RTDB_ERR) 184 if (.not.geom_rtdb_store(rtdb,geom_ref,'vscf_reference')) 185 * call errquit('vscfm: geom_rtdb_store failed',552, RTDB_ERR) 186 if (.not. rtdb_put(rtdb,'vscf:restart',mt_log,1,restart)) 187 $ call errquit('vscfm: rtdb_put restart',552, RTDB_ERR) 188 end if 189c 190c pick up a new geometry that is used for the potential scan 191c 192 if (.not.geom_create(geom_new,'geometry')) 193 * call errquit('vscfm: geom_new geom_create failed',555, GEOM_ERR) 194 if (.not.geom_rtdb_load(rtdb,geom_new,'geometry')) 195 * call errquit('vscfm: geom_rtdb_load failed',555, RTDB_ERR) 196 if (.not.geom_strip_sym(geom_new)) call errquit 197 * ('vscfm: geom_strip_sym failed',555, GEOM_ERR) 198C 199C Calculate diagonal and coupling potentials on grids 200C 201 CALL VGRID(vec,eigen,dbl_mb(k_sfreq),dbl_mb(k_rq),dbl_mb(k_dq), 202 * dbl_mb(k_dx),coords(1,1,geom_ref),dbl_mb(k_diagv), 203 * dbl_mb(k_coupv),dbl_mb(k_XTRIPV),dbl_mb(k_dmx), 204 * dbl_mb(k_dmy), 205 * dbl_mb(k_dmz),dbl_mb(k_dm2x),dbl_mb(k_dm2y), 206 * dbl_mb(k_dm2z),NC1,NAT,NNM,NGRID,NDIPMO,NSTART, 207 * NCOUP,rtdb,geom_new,ldiag,lcoup,ltrip) 208 IF (ga_nodeid().eq.0) WRITE(LuOut,9060) 209c 210c get rid of some memory that we no longer need 211c 212 if (.not. ma_pop_stack(l_dx)) 213 & call errquit('vscfm:ma_pop of l_dx failed',555, MA_ERR) 214 if (.not. ma_pop_stack(l_dm2z)) 215 & call errquit('vscfm:ma_pop of l_dm2z failed',555, MA_ERR) 216 if (.not. ma_pop_stack(l_dm2y)) 217 & call errquit('vscfm:ma_pop of l_dm2y failed',555, MA_ERR) 218 if (.not. ma_pop_stack(l_dm2x)) 219 & call errquit('vscfm:ma_pop of l_dm2x failed',555, MA_ERR) 220c 221c get memory for the next section 222c 223 if (.not. 224 & ma_push_get(mt_int,nnm,'state indices',l_state,k_state)) 225 & call errquit('vscfm: could not allocate l_state',555, MA_ERR) 226 if (.not. 227 & ma_push_get(mt_dbl,ngrid,'temp V',l_tv,k_tv)) 228 & call errquit('vscfm: could not allocate l_tv',555, MA_ERR) 229 if (.not. 230 & ma_push_get(mt_dbl,nnm*ngrid,'vscf',l_vscf,k_vscf)) 231 & call errquit('vscfm: could not allocate l_vscf',555, MA_ERR) 232 if (.not. 233 & ma_push_get(mt_dbl,ngrid,'vhf',l_vhf,k_vhf)) 234 & call errquit('vscfm: could not allocate l_vhf',555, MA_ERR) 235 if (.not. 236 & ma_push_get(mt_dbl,nnm,'enrgy',l_enrgy,k_enrgy)) 237 & call errquit('vscfm: could not allocate l_enrgy',555, MA_ERR) 238 if (.not. 239 & ma_push_get(mt_dbl,nnm+1,'ediag',l_ediag,k_ediag)) 240 & call errquit('vscfm: could not allocate l_ediag',555, MA_ERR) 241 if (.not. 242 & ma_push_get(mt_dbl,nnm+1,'escf',l_escf,k_escf)) 243 & call errquit('vscfm: could not allocate l_escf',555, MA_ERR) 244 if (.not. 245 & ma_push_get(mt_dbl,nnm+1,'emppt',l_emppt,k_emppt)) 246 & call errquit('vscfm: could not allocate l_emppt',555, MA_ERR) 247 if (.not. 248 & ma_push_get(mt_dbl,nnm*ngrid,'wave',l_wave,k_wave)) 249 & call errquit('vscfm: could not allocate l_wave',555, MA_ERR) 250 if (.not. 251 & ma_push_get(mt_dbl,nnm*ngrid*nst,'vwave',l_vwave,k_vwave)) 252 & call errquit('vscfm: could not allocate l_vwave',555, MA_ERR) 253 if (.not. 254 & ma_push_get(mt_dbl,nnm*nst,'virte',l_virte,k_virte)) 255 & call errquit('vscfm: could not allocate l_virte',555, MA_ERR) 256 if (.not. 257 & ma_push_get(mt_int,nnm*nvirst,'vst',l_vst,k_vst)) 258 & call errquit('vscfm: could not allocate l_vst',555, MA_ERR) 259 if (.not. 260 & ma_push_get(mt_int,ngrid,'pvt',l_pvt,k_pvt)) 261 & call errquit('vscfm: could not allocate l_pvt',555, MA_ERR) 262 if (.not. 263 & ma_push_get(mt_dbl,ngrid,'xx',l_xx,k_xx)) 264 & call errquit('vscfm: could not allocate l_xx',555, MA_ERR) 265 if (.not. 266 & ma_push_get(mt_dbl,ngrid,'a',l_a,k_a)) 267 & call errquit('vscfm: could not allocate l_a',555, MA_ERR) 268 if (.not. 269 & ma_push_get(mt_dbl,ngrid3,'phi',l_phi,k_phi)) 270 & call errquit('vscfm: could not allocate l_phi',555, MA_ERR) 271 if (.not. 272 & ma_push_get(mt_dbl,ngrid3,'r',l_r,k_r)) 273 & call errquit('vscfm: could not allocate l_r',555, MA_ERR) 274 if (.not. 275 & ma_push_get(mt_dbl,ngrid3,'rr',l_rr,k_rr)) 276 & call errquit('vscfm: could not allocate l_rr',555, MA_ERR) 277 if (.not. 278 & ma_push_get(mt_dbl,ngrid3,'g',l_g,k_g)) 279 & call errquit('vscfm: could not allocate l_g',555, MA_ERR) 280 if (.not. 281 & ma_push_get(mt_dbl,ngrid3,'v',l_v,k_v)) 282 & call errquit('vscfm: could not allocate l_v',555, MA_ERR) 283 if (.not. 284 & ma_push_get(mt_dbl,ngrid2,'h',l_h,k_h)) 285 & call errquit('vscfm: could not allocate l_h',555, MA_ERR) 286 if (.not. 287 & ma_push_get(mt_dbl,ngrid,'ec',l_ec,k_ec)) 288 & call errquit('vscfm: could not allocate l_ec',555, MA_ERR) 289 if (.not. 290 & ma_push_get(mt_dbl,ngrid*ngrid,'vecc',l_vecc,k_vecc)) 291 & call errquit('vscfm: could not allocate l_vecc',555, MA_ERR) 292 if (.not. 293 & ma_push_get(mt_dbl,ngrid*8,'scr1',l_scr1,k_scr1)) 294 & call errquit('vscfm: could not allocate l_scr1',555, MA_ERR) 295 if (.not. 296 & ma_push_get(mt_int,ngrid,'ia1',l_ia1,k_ia1)) 297 & call errquit('vscfm: could not allocate l_ia1',555, MA_ERR) 298 if (.not. 299 & ma_push_get(mt_dbl,ngrid3,'gr',l_gr,k_gr)) 300 & call errquit('vscfm: could not allocate l_gr',555, MA_ERR) 301 if (.not. 302 & ma_push_get(mt_dbl,ngrid,'tmp',l_tmp,k_tmp)) 303 & call errquit('vscfm: could not allocate l_tmp',555, MA_ERR) 304 if (.not. 305 & ma_push_get(mt_dbl,ngrid,'twave',l_twave,k_twave)) 306 & call errquit('vscfm: could not allocate l_twave',555, MA_ERR) 307 if (.not. 308 & ma_push_get(mt_dbl,nvirst,'emp0',l_emp0,k_emp0)) 309 & call errquit('vscfm: could not allocate l_emp0',555, MA_ERR) 310 if (.not. 311 & ma_push_get(mt_dbl,nvirst,'vmp',l_vmp,k_vmp)) 312 & call errquit('vscfm: could not allocate l_vmp',555, MA_ERR) 313 if (.not. 314 & ma_push_get(mt_dbl,nnm,'ovrlp',l_ovrlp,k_ovrlp)) 315 & call errquit('vscfm: could not allocate l_ovrlp',555, MA_ERR) 316 if (.not. 317 & ma_push_get(mt_int,nvirst*nnm,'virt',l_virt,k_virt)) 318 & call errquit('vscfm: could not allocate l_virt',555, MA_ERR) 319 if (.not. 320 & ma_push_get(mt_int,nnm,'ref',l_ref,k_ref)) 321 & call errquit('vscfm: could not allocate l_ref',555, MA_ERR) 322C 323C Perform VSCF and MP2-VSCF 324C 325 IF (ga_nodeid().eq.0) THEN 326 CALL VSCFMP(dbl_mb(k_sfreq),dbl_mb(k_rq),dbl_mb(k_dq), 327 * dbl_mb(k_diagv),dbl_mb(k_coupv),int_mb(k_state), 328 * dbl_mb(k_tv),dbl_mb(k_vscf),dbl_mb(k_vhf), 329 * dbl_mb(k_enrgy),dbl_mb(k_ediag),dbl_mb(k_escf), 330 * dbl_mb(k_emppt),dbl_mb(k_wave),dbl_mb(k_awave), 331 * dbl_mb(k_vwave),dbl_mb(k_virte),int_mb(k_vst), 332 * int_mb(k_pvt),dbl_mb(k_xx),dbl_mb(k_a), 333 * dbl_mb(k_phi),dbl_mb(k_r),dbl_mb(k_rr), 334 * dbl_mb(k_g),dbl_mb(k_v),dbl_mb(k_h), 335 * dbl_mb(k_ec),dbl_mb(k_vecc),dbl_mb(k_scr1), 336 * int_mb(k_ia1),dbl_mb(k_gr),dbl_mb(k_tmp), 337 * dbl_mb(k_twave),dbl_mb(k_emp0),dbl_mb(k_vmp), 338 * dbl_mb(k_ovrlp),int_mb(k_virt),int_mb(k_ref), 339 * dbl_mb(k_fvscf),dbl_mb(k_fmp2),dbl_mb(k_XTRIPV), 340 * vcfct,NNM,NGRID,NGRID2,NST,NVIRST,ncoup,nmax,iexc) 341 WRITE(LuOut,9070) 342 END IF 343c 344c get rid of some more memory that isn't needed 345c 346 if (.not. ma_pop_stack(l_ref)) 347 & call errquit('vscfm:ma_pop of l_ref failed',555, MA_ERR) 348 if (.not. ma_pop_stack(l_virt)) 349 & call errquit('vscfm:ma_pop of l_virt failed',555, MA_ERR) 350 if (.not. ma_pop_stack(l_ovrlp)) 351 & call errquit('vscfm:ma_pop of l_ovrlp failed',555, MA_ERR) 352 if (.not. ma_pop_stack(l_vmp)) 353 & call errquit('vscfm:ma_pop of l_vmp failed',555, MA_ERR) 354 if (.not. ma_pop_stack(l_emp0)) 355 & call errquit('vscfm:ma_pop of l_emp0 failed',555, MA_ERR) 356 if (.not. ma_pop_stack(l_twave)) 357 & call errquit('vscfm:ma_pop of l_twave failed',555, MA_ERR) 358 if (.not. ma_pop_stack(l_tmp)) 359 & call errquit('vscfm:ma_pop of l_tmp failed',555, MA_ERR) 360 if (.not. ma_pop_stack(l_gr)) 361 & call errquit('vscfm:ma_pop of l_gr failed',555, MA_ERR) 362 if (.not. ma_pop_stack(l_ia1)) 363 & call errquit('vscfm:ma_pop of l_ia1 failed',555, MA_ERR) 364 if (.not. ma_pop_stack(l_scr1)) 365 & call errquit('vscfm:ma_pop of l_scr1 failed',555, MA_ERR) 366 if (.not. ma_pop_stack(l_vecc)) 367 & call errquit('vscfm:ma_pop of l_vecc failed',555, MA_ERR) 368 if (.not. ma_pop_stack(l_ec)) 369 & call errquit('vscfm:ma_pop of l_ec failed',555, MA_ERR) 370 if (.not. ma_pop_stack(l_h)) 371 & call errquit('vscfm:ma_pop of l_h failed',555, MA_ERR) 372 if (.not. ma_pop_stack(l_v)) 373 & call errquit('vscfm:ma_pop of l_v failed',555, MA_ERR) 374 if (.not. ma_pop_stack(l_g)) 375 & call errquit('vscfm:ma_pop of l_g failed',555, MA_ERR) 376 if (.not. ma_pop_stack(l_rr)) 377 & call errquit('vscfm:ma_pop of l_rr failed',555, MA_ERR) 378 if (.not. ma_pop_stack(l_r)) 379 & call errquit('vscfm:ma_pop of l_r failed',555, MA_ERR) 380 if (.not. ma_pop_stack(l_phi)) 381 & call errquit('vscfm:ma_pop of l_phi failed',555, MA_ERR) 382 if (.not. ma_pop_stack(l_a)) 383 & call errquit('vscfm:ma_pop of l_a failed',555, MA_ERR) 384 if (.not. ma_pop_stack(l_xx)) 385 & call errquit('vscfm:ma_pop of l_xx failed',555, MA_ERR) 386 if (.not. ma_pop_stack(l_pvt)) 387 & call errquit('vscfm:ma_pop of l_pvt failed',555, MA_ERR) 388 if (.not. ma_pop_stack(l_vst)) 389 & call errquit('vscfm:ma_pop of l_vst failed',555, MA_ERR) 390 if (.not. ma_pop_stack(l_virte)) 391 & call errquit('vscfm:ma_pop of l_virte failed',555, MA_ERR) 392 if (.not. ma_pop_stack(l_vwave)) 393 & call errquit('vscfm:ma_pop of l_vwave failed',555, MA_ERR) 394 if (.not. ma_pop_stack(l_wave)) 395 & call errquit('vscfm:ma_pop of l_wave failed',555, MA_ERR) 396 if (.not. ma_pop_stack(l_emppt)) 397 & call errquit('vscfm:ma_pop of l_emppt failed',555, MA_ERR) 398 if (.not. ma_pop_stack(l_escf)) 399 & call errquit('vscfm:ma_pop of l_escf failed',555, MA_ERR) 400 if (.not. ma_pop_stack(l_ediag)) 401 & call errquit('vscfm:ma_pop of l_ediag failed',555, MA_ERR) 402 if (.not. ma_pop_stack(l_enrgy)) 403 & call errquit('vscfm:ma_pop of l_enrgy failed',555, MA_ERR) 404 if (.not. ma_pop_stack(l_vhf)) 405 & call errquit('vscfm:ma_pop of l_vhf failed',555, MA_ERR) 406 if (.not. ma_pop_stack(l_vscf)) 407 & call errquit('vscfm:ma_pop of l_vscf failed',555, MA_ERR) 408 if (.not. ma_pop_stack(l_tv)) 409 & call errquit('vscfm:ma_pop of l_tv failed',555, MA_ERR) 410 if (.not. ma_pop_stack(l_state)) 411 & call errquit('vscfm:ma_pop of l_state failed',555, MA_ERR) 412 if (.not. ma_pop_stack(l_coupv)) 413 & call errquit('vscfm:ma_pop of l_coupv failed',555, MA_ERR) 414 if (.not. ma_pop_stack(l_diagv)) 415 & call errquit('vscfm:ma_pop of l_diagv failed',555, MA_ERR) 416 if (.not. ma_pop_stack(l_sfreq)) 417 & call errquit('vscfm:ma_pop of l_sfreq failed',555, MA_ERR) 418c 419c get memory for the next section 420c 421 if (.not. 422 & ma_push_get(mt_dbl,nnm,'intensity',l_int,k_int)) 423 & call errquit('vscfm: could not allocate l_int',555, MA_ERR) 424C 425C Calculate IR intensities 426C 427 if (ga_nodeid().eq.0) then 428 IF (ndipmo.eq.0) THEN 429 if (.not. 430 & ma_push_get(mt_dbl,nnm,'dder',l_dder,k_dder)) 431 & call errquit('vscfm: could not allocate l_dder',555, MA_ERR) 432 CALL DINTENS(dbl_mb(k_int),dbl_mb(k_ddm),dbl_mb(k_dder), 433 * dbl_mb(k_fvscf),dbl_mb(k_fmp2),vec,dbl_mb(k_rq), 434 * dbl_mb(k_dq),dbl_mb(k_awave),NNM,NGRID,NC1) 435 if (.not. ma_pop_stack(l_dder)) 436 & call errquit('vscfm:ma_pop of l_dder failed',555, MA_ERR) 437 ELSE 438 CALL INTENS(dbl_mb(k_int),dbl_mb(k_dmx),dbl_mb(k_dmy), 439 * dbl_mb(k_dmz),dbl_mb(k_fvscf),dbl_mb(k_fmp2), 440 * dbl_mb(k_dq),dbl_mb(k_awave),NNM,NGRID) 441 END IF 442 endif 443c 444c clean up memory 445c 446 if (.not. ma_pop_stack(l_int)) 447 & call errquit('vscfm:ma_pop of l_int failed',555, MA_ERR) 448 if (ndipmo.eq.0) then 449 if (.not. ma_pop_stack(l_ddm)) 450 & call errquit('vscfm:ma_pop of l_dmz failed',555, MA_ERR) 451 else 452 if (.not. ma_pop_stack(l_dmz)) 453 & call errquit('vscfm:ma_pop of l_dmz failed',555, MA_ERR) 454 if (.not. ma_pop_stack(l_dmy)) 455 & call errquit('vscfm:ma_pop of l_dmy failed',555, MA_ERR) 456 if (.not. ma_pop_stack(l_dmx)) 457 & call errquit('vscfm:ma_pop of l_dmx failed',555, MA_ERR) 458 endif 459 if (.not. ma_pop_stack(l_awave)) 460 & call errquit('vscfm:ma_pop of l_awave failed',555, MA_ERR) 461 if (.not. ma_pop_stack(l_rq)) 462 & call errquit('vscfm:ma_pop of l_rq failed',555, MA_ERR) 463 if (.not. ma_pop_stack(l_dq)) 464 & call errquit('vscfm:ma_pop of l_dq failed',555, MA_ERR) 465 if (.not. ma_pop_stack(l_fmp2)) 466 & call errquit('vscfm:ma_pop of l_fmp2 failed',555, MA_ERR) 467 if (.not. ma_pop_stack(l_fvscf)) 468 & call errquit('vscfm:ma_pop of l_fvscf failed',555, MA_ERR) 469 if (.not. ma_pop_stack(l_XTRIPV)) 470 & call errquit('vscfm:ma_pop of l_XTRIPV failed',555, MA_ERR) 471c 472c Reset geometry to default (reference) geometry 473c 474 if (.not.geom_rtdb_load(rtdb,geom_ref,'vscf_reference')) 475 * call errquit('vscfm: geom_rtdb_load failed',221, RTDB_ERR) 476 if (.not.geom_rtdb_store(rtdb,geom_ref,'geometry')) 477 * call errquit('vscfm: geom_rtdb_store failed',221, RTDB_ERR) 478 if (.not.geom_rtdb_delete(rtdb,'vscf_reference')) 479 * call errquit('vscfm: geom_rtdb_store failed',221, RTDB_ERR) 480 if (.not.geom_destroy(geom_new)) 481 * call errquit('vscfm: geom_destroy',221, RTDB_ERR) 482 if (.not.geom_destroy(geom_ref)) 483 * call errquit('vscfm: geom_destroy',221, RTDB_ERR) 484c 485 IF (ga_nodeid().eq.0) WRITE(LuOut,9080) 486 return 487c 48889900 call errquit('vscfm: cannot open fd_ddipole file',221, DISK_ERR) 48989901 call errquit('vscfm: error reading fd_ddipole file',221, DISK_ERR) 49089902 call errquit('vscfm: incomplete fd_ddipole file',221, DISK_ERR) 491 9060 FORMAT(1X,'......Done with potentials on grids......') 492 9070 FORMAT(1X,'......Finished vibrational SCF......') 493 9080 FORMAT(1X,'......Finished calculating IR intensities......') 494 end 495c*module vscf *deck vgrid 496 SUBROUTINE VGRID(VEC,EIG,FREQ,RQ,DQ,DX,C0,DIAGV,COUPV,TRIPV, 497 * DM1X,DM1Y,DM1Z,DM2X,DM2Y,DM2Z, 498 * NC1,NAT1,NNM,NGRID,NDIPMO, 499 * NFIN,NCOUP,rtdb,geom, 500 * ldiag,lcoup,ltrip) 501C 502 implicit none 503#include "errquit.fh" 504c 505#include "nwc_const.fh" 506#include "geomP.fh" 507#include "geom.fh" 508#include "global.fh" 509#include "mafdecls.fh" 510#include "rtdb.fh" 511#include "util.fh" 512#include "stdio.fh" 513C 514 integer NC1,NAT1,NNM,NGRID,NDIPMO,NCOUP 515 LOGICAL psi0,status 516 integer ldiag, lcoup, ltrip ! these are unit numbers 517 character*32 theory 518 character*255 vectors_in 519c 520 logical task_energy, property 521 external task_energy, property 522C 523 integer rtdb, geom 524c 525 double precision EIG(NC1),VEC(NC1,NC1),DX(NC1),DIAGV(NNM,NGRID) 526 double precision FREQ(NNM),C0(3,NAT1),COUPV(NNM,NNM,NGRID,NGRID) 527 double precision DQ(NNM),RQ(NNM,NGRID),TRIPV(NNM,NNM,NNM, 528 * NGRID,NGRID,NGRID) 529 double precision DM1X(NNM,NGRID),DM1Y(NNM,NGRID),DM1Z(NNM,NGRID) 530 double precision DM2X(NNM,NNM,NGRID,NGRID), 531 * DM2Y(NNM,NNM,NGRID,NGRID), 532 * DM2Z(NNM,NNM,NGRID,NGRID) 533 double precision dmtemp(3), dmx, dmy, dmz, dmx0, dmy0, dmz0 534 double precision qmin, qrange, wave, fact, e, e0 535 integer NFIN(7),nstart, n, m, j, i, k 536 integer mn, jm, km, jl, kl, im, il, jk, ii 537C 538 double precision zero, one, four, two, amu 539 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, FOUR=4.0D+00) 540 PARAMETER (TWO=2.0D+00, AMU=1.8229D+03) 541C 542 FACT = one/sqrt(AMU) 543 NSTART=NC1-NNM+1 544c 545c mp2 runs must execute gradient code to get density matrix 546c unfortunately, parallel runs don't have fast cutout yet. 547c 548c if(mplevl.gt.0) then 549c mpprop=1 550c if(nproc.gt.1) then 551c if(ga_nodeid().eq.0) 552c & write(LuOut,*) 'need parallel property cutout' 553c call errquit('vgrid: need mp2 nos for mp2 properties',555) 554c end if 555c end if 556C 557c do we need to get starting orbitals at the input geometry? 558c 559 psi0 = .true. 560c 561 if(.not.rtdb_cget(rtdb,'task:theory',1,theory)) 562 + call errquit('vgrid: no input for theory?',0, RTDB_ERR) 563 if(psi0) then 564 if(ga_nodeid().eq.0) write(LuOut,9005) 565 if (.not.task_energy(rtdb)) ! equilibrium geometry energy 566 & call errquit('vgrid: first energy did not converge!',555, 567 & CALC_ERR) 568 if (theory.eq.'dft') then 569 if (.not.rtdb_get(rtdb,'dft:energy',mt_dbl,1,e)) 570 & call errquit('vgrid: failed to get energy from rtdb 1',555, 571 & RTDB_ERR) 572 else 573 if (.not.rtdb_get(rtdb,'scf:energy',mt_dbl,1,e)) 574 & call errquit('vgrid: failed to get energy from rtdb 1',555, 575 & RTDB_ERR) 576 endif 577c 578c Point prop:input vectors to the scf or dft vectors so it does 579c not have to redo the scf/dft everytime 580c 581 if (theory.eq.'dft') then 582 if (.not. rtdb_cget(rtdb,'dft:output vectors',1,vectors_in)) 583 & call errquit('vgrid: rtdb_cget dft:output failed', 100, 584 & RTDB_ERR) 585 else 586 if (.not. rtdb_cget(rtdb,'scf:output vectors',1,vectors_in)) 587 & call errquit('vgrid: rtdb_cget scf:output failed', 100, 588 & RTDB_ERR) 589 end if 590 if (.not. rtdb_cput(rtdb,'prop:vectors',1,vectors_in)) 591 & call errquit('vgrid: rtdb_cput prop:input failed', 100, 592 & RTDB_ERR) 593c 594 if (.not.rtdb_put(rtdb,'prop:dipole',mt_int,1,0)) 595 & call errquit('vgrid: rtdb_put of dipole failed',555, 596 & RTDB_ERR) 597 if (.not.property(rtdb)) ! equilibrium dipole moment 598 & call errquit('vgrid: first dipole did not work!',555, 599 & RTDB_ERR) 600 if (.not.rtdb_get(rtdb,'prop:dipval',mt_dbl,3,dmtemp)) 601 & call errquit('vgrid: rtdb_get of dipval failed',555, 602 & RTDB_ERR) 603 dmx = dmtemp(1) 604 dmy = dmtemp(2) 605 dmz = dmtemp(3) 606 else 607 e = zero 608 dmx = zero 609 dmy = zero 610 dmz = zero 611 end if 612C 613C save potential energy, dipole moment at equilibrium 614C 615 E0=E 616 dmx0=dmx 617 dmy0=dmy 618 dmz0=dmz 619C 620C Frequencies in atomic units 621C 622 WAVE = 1.6021892D00*1.6021892D00*1.0D2 623 WAVE = WAVE/(8.0d0*ACOS(0.0D00)*8.85418782D00*0.5291771D00) 624 WAVE = WAVE*6.022045D00/(0.5291771D00*0.5291771D00) 625 WAVE = SQRT(WAVE) 626 WAVE = WAVE/(4.0d0*ACOS(0.0D00)*2.9979245D00) 627 WAVE = WAVE*1.0D4 628c 629c Bring in frequencies and convert back from cm-1, and then convert to atomic units 630c 631 do i=NC1,NSTART,-1 632 ii=NC1+1-i 633 freq(ii)=(eig(i)/WAVE)/sqrt(AMU) 634 end do 635c 636c If all diag values were done, goto coupl 637c 638 if (nfin(1).gt.0) goto 100 639C 640C Calculate diagonal V on a grid 641C 642 if (ga_nodeid().eq.0) write(LuOut,9040) ngrid, nc1-nstart+1 643C 644 do i=NC1,NSTART,-1 645 im=nc1-i+1 646 Qrange=FOUR 647 Qrange=Qrange/sqrt(freq(im)) 648 Qmin=-Qrange 649 dQ(im)=(TWO*Qrange)/(ngrid-1) 650 do il=1, ngrid 651 RQ(im,il)=Qmin+(il-1)*dQ(im) 652 jk= 0 653 do j=1, nat1 654 do k=1, 3 655 jk=jk+1 656 dx(jk) = FACT*VEC(jk,i)*RQ(im,il) 657 coords(k,j,geom) = C0(k,j) + dx(jk) 658 end do 659 end do 660c 661c Skip pieces on restart until restart point is reached 662c At that point restart point reset to zero to indicate we 663c just want to continue to calculate points 664c 665c If nfin(2).gt.zero, then we have not reached the restart 666c point and we need to skip one more 667c 668 if (im.eq.nfin(2).and.il.eq.nfin(3)) then 669 nfin(2)=0 670 goto 99 671 else if (nfin(2).gt.0) then 672 goto 99 673 endif 674c 675 if (.not.geom_rtdb_store(rtdb,geom,'geometry')) 676 * call errquit('vgrid: geom_rtdb_store failed',555, 677 & RTDB_ERR) 678 if(ga_nodeid().eq.0) write(LuOut,9045) il,i 679 status = task_energy(rtdb) 680 if (theory.eq.'dft') then 681 if (.not.rtdb_get(rtdb,'dft:energy',mt_dbl,1,e)) 682 & call errquit('vgrid: failed to get energy from rtdb 2',555, 683 & RTDB_ERR) 684 else 685 if (.not.rtdb_get(rtdb,'scf:energy',mt_dbl,1,e)) 686 & call errquit('vgrid: failed to get energy from rtdb 2',555, 687 & RTDB_ERR) 688 endif 689 IF(.not.status .AND. ga_nodeid().eq.0) WRITE(LuOut,9000) 690c 691 diagV(im,il) = e - e0 692 if (ndipmo.eq.1) then 693 if (status) then 694 if (.not.property(rtdb)) 695 & call errquit('vgrid: first dipole did not work!',555, 696 & RTDB_ERR) 697 if (.not.rtdb_get(rtdb,'prop:dipval',mt_dbl,3,dmtemp)) 698 & call errquit('vgrid: rtdb_get of dipval failed',555, 699 & RTDB_ERR) 700 dmx = dmtemp(1) 701 dmy = dmtemp(2) 702 dmz = dmtemp(3) 703 else 704 dmx = dmx0 705 dmy = dmy0 706 dmz = dmz0 707 end if 708 dm1x(im,il) = dmx - dmx0 709 dm1y(im,il) = dmy - dmy0 710 dm1z(im,il) = dmz - dmz0 711 if (ga_nodeid().eq.0) write(ldiag,9015) im, RQ(im,il), 712 * diagV(im,il),dm1x(im,il),dm1y(im,il),dm1z(im,il) 713 if(ga_nodeid().eq.0) call util_flush(ldiag) 714 else 715 if (ga_nodeid().eq.0) 716 & write(ldiag,9010) im, RQ(im,il), diagV(im,il) 717 if(ga_nodeid().eq.0) call util_flush(ldiag) 718 end if 719 99 continue 720 end do 721 end do 722 if (ga_nodeid().eq.0) write(LuOut,9050) 723C 724C Off-diagonal pair coupling potential 725C 726 100 CONTINUE 727 if ((ncoup.le.1) .or. (nfin(1).gt.1)) goto 200 728C 729 if (ga_nodeid().eq.0) write(LuOut,9060) ngrid,ngrid, 730 * (nc1-nstart+1)*(nc1-nstart)/2 731C 732 do i=NC1, NSTART+1,-1 733 im=nc1-i+1 734 do j=i-1, NSTART,-1 735 jm=nc1-j+1 736 do il=1, ngrid 737 do jl=1, ngrid 738 mn=0 739 do m=1, nat1 740 do n=1, 3 741 mn=mn+1 742 dx(mn)= FACT*VEC(mn,i)*RQ(im,il) + 743 * FACT*VEC(mn,j)*RQ(jm,jl) 744 coords(n,m,geom) = C0(n,m) + dx(mn) 745 end do 746 end do 747c 748c Skip pieces on restart until restart point is reached 749c At that point restart point reset to zero to indicate we 750c just want to continue to calculate points 751c 752c If nfin(2).gt.zero, then we have not reached the restart 753c point and we need to skip one more 754c 755 if (nfin(1) .le. 1) then ! Do them all 756 else if (im.eq.nfin(2).and.jm.eq.nfin(4).and. 757 & il.eq.nfin(3).and.jl.eq.nfin(5)) then 758 nfin(2)=0 759 goto 199 760 else if (nfin(2).gt.0) then 761 goto 199 762 endif 763c 764 if (.not.geom_rtdb_store(rtdb,geom,'geometry')) 765 * call errquit('vgrid: geom_rtdb_store failed',555, 766 & RTDB_ERR) 767 if(ga_nodeid().eq.0) write(LuOut,9065) il,jl,i,j 768 status = task_energy(rtdb) 769 if (theory.eq.'dft') then 770 if (.not.rtdb_get(rtdb,'dft:energy',mt_dbl,1,e)) 771 * call errquit 772 * ('vgrid: failed to get energy from rtdb 3',555, 773 & RTDB_ERR) 774 else 775 if (.not.rtdb_get(rtdb,'scf:energy',mt_dbl,1,e)) 776 * call errquit 777 * ('vgrid: failed to get energy from rtdb 3',555, 778 & RTDB_ERR) 779 endif 780 IF(.not.status .AND. ga_nodeid().eq.0) 781 * WRITE(LuOut,9000) 782c 783 coupV(im,jm,il,jl)=e-diagV(im,il)-diagV(jm,jl)-e0 784 if (ndipmo.eq.1) then 785 if (status) then 786 if (.not.property(rtdb)) 787 & call errquit('vgrid: dipole did not work!',555, 788 & RTDB_ERR) 789c 790c Need to not errquit in above line and write an error message 791c 792 if (.not.rtdb_get(rtdb,'prop:dipval', 793 & mt_dbl,3,dmtemp)) 794 & call errquit 795 & ('vgrid: rtdb_get of dipval failed',555, RTDB_ERR) 796 dmx = dmtemp(1) 797 dmy = dmtemp(2) 798 dmz = dmtemp(3) 799 else 800 dmx = dmx0 801 dmy = dmy0 802 dmz = dmz0 803 end if 804 dm2x(im,jm,il,jl) = dmx - dmx0 805 dm2y(im,jm,il,jl) = dmy - dmy0 806 dm2z(im,jm,il,jl) = dmz - dmz0 807 if (ga_nodeid().eq.0) then 808#if defined(FUJITSU_SOLARIS) || defined(SOLARIS) || defined(GCC46) 809 backspace 82 810#endif 811 write(lcoup,9022) im, jm, RQ(im,il), RQ(jm,jl), 812 * coupV(im,jm,il,jl),dm2x(im,jm,il,jl), 813 * dm2y(im,jm,il,jl),dm2z(im,jm,il,jl) 814 call util_flush(lcoup) 815 endif 816 else 817 if(ga_nodeid().eq.0) 818 & write(lcoup,9020) im,jm,RQ(im,il),RQ(jm,jl), 819 * coupV(im,jm,il,jl) 820 if(ga_nodeid().eq.0) call util_flush(lcoup) 821 end if 822 199 continue 823 end do 824 end do 825 end do 826 end do 827c 828 e=e0 829 if (ga_nodeid().eq.0) write(LuOut,9070) 830C 831C 3-body coupling potential 832C 833 200 CONTINUE 834 if ((ncoup.le.2) .or. (nfin(1).eq.3)) goto 300 835C 836 if (ga_nodeid().eq.0) then 837 write(LuOut,*) 838 write(LuOut,*) ' 3-body coupling potential' 839 write(LuOut,*) 840 end if 841C 842 if (ga_nodeid().eq.0) write(LuOut,9080) 843 844 do i=NC1, NSTART+2,-1 845 im=nc1-i+1 846 do j=i-1, NSTART+1,-1 847 jm=nc1-j+1 848 do k=j-1, NSTART,-1 849 km=nc1-k+1 850 do il=1, ngrid 851 do jl=1, ngrid 852 do kl=1, ngrid 853 mn=0 854 do m=1, nat1 855 do n=1, 3 856 mn=mn+1 857 dx(mn)= FACT*VEC(mn,i)*RQ(im,il) + 858 * FACT*VEC(mn,j)*RQ(jm,jl) + 859 * FACT*VEC(mn,k)*RQ(km,kl) 860 coords(n,m,geom) = C0(n,m) + dx(mn) 861 end do 862 end do 863c 864c Skip pieces on restart until restart point is reached 865c At that point restart point reset to zero to indicate we 866c just want to continue to calculate points 867c 868c If nfin(2).gt.zero, then we have not reached the restart 869c point and we need to skip one more 870c 871 if (nfin(1) .le. 2) then ! Do them all 872 else if (im.eq.nfin(2).and.jm.eq.nfin(4).and. 873 & il.eq.nfin(3).and.jl.eq.nfin(5).and. 874 & km.eq.nfin(6).and.kl.eq.nfin(7)) then 875 nfin(2)=0 876 goto 299 877 else if (nfin(2).gt.0) then 878 goto 299 879 endif 880c 881 if (.not.geom_rtdb_store(rtdb,geom,'geometry')) 882 * call errquit('vgrid: geom_rtdb_store failed',555, 883 * RTDB_ERR) 884 if(ga_nodeid().eq.0) 885 * write(LuOut,9085) il,jl,kl,i,j,k 886 status = task_energy(rtdb) 887 if (theory.eq.'dft') then 888 if (.not.rtdb_get(rtdb,'dft:energy',mt_dbl,1,e)) 889 * call errquit 890 * ('vgrid: failed to get energy from rtdb',555, 891 * RTDB_ERR) 892 else 893 if (.not.rtdb_get(rtdb,'scf:energy',mt_dbl,1,e)) 894 * call errquit 895 * ('vgrid: failed to get energy from rtdb',555, 896 * RTDB_ERR) 897 endif 898 IF(.not.status .AND. ga_nodeid().eq.0) 899 * WRITE(LuOut,9000) 900 901 tripV(im,jm,km,il,jl,kl)=e-e0- 902 * diagV(im,il)-diagV(jm,jl)-diagV(km,kl)- 903 * coupV(im,jm,il,jl)-coupV(im,km,il,kl)- 904 * coupV(jm,km,jl,kl) 905 if(ga_nodeid().eq.0) write(LuOut,9025) im,jm,km, 906 * RQ(im,il),RQ(jm,jl),RQ(km,kl), 907 * tripV(im,jm,km,il,jl,kl) 908 if(ga_nodeid().eq.0) call util_flush(LuOut) 909 299 continue 910 end do 911 end do 912 end do 913 end do 914 end do 915 end do 916 917 e=e0 918 if (ga_nodeid().eq.0) write(LuOut,9090) 919c 920 300 CONTINUE 921 if (.not. rtdb_delete(rtdb,'prop:vectors')) 922 & call errquit('vgrid: rtdb_delete failed',100,RTDB_ERR) 923c 924 RETURN 925C 926 9000 FORMAT(1x,4(2h*-),'*'/1X,'warning !'/1x,4(2h*-),'*'/ 927 * 1x,'SCF HAS NOT CONVERGED at this vscf point!') 928 9005 format(//1x,'evaluating wavefunction at original geometry...'//) 929 9010 format(1x,I2,2x,f14.8,2x,e14.7) 930 9015 format(1x,I2,2x,f14.8,2x,e14.7,3(2x,e14.7)) 931 9020 format(1x,2(I2,2x),2(f14.8,2x),e14.7) 932 9022 format(1x,2(I2,2x),2(f14.8,2x),e14.7,3(e14.7,2x)) 933 9025 format(1x,3(I2,2x),3(f14.8,2x),e14.7) 934 9030 format(1x,f14.9) 935 9040 FORMAT(//1X,'Starting diagonal potential on a grid of ',I3, 936 * ' points for ',I3,' normal modes') 937 9045 format(//1x,'VSCF: Energy and dipole for grid point',i3, 938 * ' along mode',i3) 939 9050 FORMAT(/1X,'Done with diagonal potential') 940 9060 FORMAT(//1X,'Starting pair coupling potential on a square grid'/ 941 * ' of',I3,' by',I3,' points for',I6,' pairs of normal modes') 942 9065 format(//1x,'VSCF: Energy and dipole for pair grid points',2i3, 943 * ' for mode pair',2i3) 944 9070 FORMAT(/1X,'Done with pair coupling potential') 945 9080 FORMAT(//1X,'Starting 3-body coupling potential') 946 9085 format(//1x,'VSCF: Energy and dipole for 3-body grid points',3i3, 947 * ' for mode triplet',3i3) 948 9090 FORMAT(/1X,'Done with 3-body coupling potential') 949 END 950C 951 SUBROUTINE VSCFMP(FREQ,RQ,DQ,DIAGV,COUPV,ISTATE,TV,VSCF,VHF, 952 * E,EDIAG,ESCF,EMPPT,WAVE,ALLWAVE,VIRTWAVE,VIRTE, 953 * IVST,IPVT,XX,A,PHI,R,RR,G,V,H,EC,VEC,SCR,IA,GR, 954 * TEMP,TWAVE,EMP0,VMP,OVRLP,JVIRT,JREF,frscf,frmp2, 955 * tripv,vcfct,NNM,NGRID,NGRID2,NST,NVST,ncoup,nmax, 956 * iexc) 957C 958 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 959#include "errquit.fh" 960c 961#include "global.fh" 962#include "stdio.fh" 963C 964 LOGICAL GOPARR,DSKWRK,MASWRK 965C 966 DIMENSION DIAGV(NNM,NGRID),COUPV(NNM,NNM,NGRID,NGRID) 967 DIMENSION TRIPV(NNM,NNM,NNM,NGRID,NGRID,NGRID) 968 DIMENSION FREQ(NNM),DQ(NNM),RQ(NNM,NGRID) 969 DIMENSION tV(NGRID),vscf(NNM,NGRID),VHF(NGRID),E(NNM) 970 dimension ediag(NNM+1), escf(NNM+1), emppt(NNM+1) 971 dimension wave(NNM,NGRID),allwave(NNM+1,NNM,NGRID) 972 dimension virtwave(NNM,NST,NGRID),virtE(NNM,NST) 973 dimension ISTATE(NNM),IVST(NVST,NNM),ipvt(NGRID) 974 dimension xx(NGRID),A(NGRID),phi(NGRID,NGRID) 975 dimension R(NGRID,NGRID),RR(NGRID,NGRID),G(NGRID,NGRID) 976 dimension V(NGRID,NGRID),H(NGRID2),EC(NGRID),VEC(NGRID,NGRID) 977 dimension GR(NGRID,NGRID),temp(NGRID),SCR(NGRID,8),IA(NGRID) 978 dimension twave(NGRID),emp0(NVST),vmp(NVST) 979 dimension jvirt(NVST,NNM),jref(NNM),ovrlp(NNM) 980 dimension frscf(nnm),frmp2(nnm) 981C 982 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) 983 PARAMETER (EPS=1.0D-06, cm=2.19474D+05) 984 PARAMETER (mxiter=100) 985 logical print_wv 986c 987c Print vibrational wave functions 988c 989 print_wv=.false. 990C 991C SCALE PAIR COUPLING POTENTIAL IF NEEDED 992C 993 vcfct1=one 994 500 fct=vcfct/vcfct1 995 if (ncoup.gt.1) then 996 do i=1, NNM-1 997 do j=i+1, NNM 998 do l=1, NGRID 999 do m=1, NGRID 1000C 1001 if ( (abs(coupV(i,j,l,m)).ge.diagV(i,l)).or. 1002 * (abs(coupV(i,j,l,m)).ge.diagV(j,m)) ) 1003 * coupV(i,j,l,m) = fct*coupV(i,j,l,m) 1004C 1005 coupV(j,i,m,l) = coupV(i,j,l,m) 1006C 1007 end do 1008 end do 1009 end do 1010 end do 1011 end if 1012C 1013C SCALE 3-BODY COUPLING POTENTIAL IF NEEDED 1014C 1015 if (ncoup.gt.2) then 1016 do i=1, NNM-2 1017 do j=i+1, NNM-1 1018 do k=j+1, NNM 1019 do l1=1, NGRID 1020 do l2=1, NGRID 1021 do l3=1, NGRID 1022C 1023 if ( (abs(tripV(i,j,k,l1,l2,l3)).ge.diagV(i,l1)).or. 1024 * (abs(tripV(i,j,k,l1,l2,l3)).ge.diagV(j,l2)).or. 1025 * (abs(tripV(i,j,k,l1,l2,l3)).ge.diagV(k,l3)) ) 1026 * tripV(i,j,k,l1,l2,l3) = fct*tripV(i,j,k,l1,l2,l3) 1027C 1028 tripV(i,k,j,l1,l3,l2) = tripV(i,j,k,l1,l2,l3) 1029 tripV(j,i,k,l2,l1,l3) = tripV(i,j,k,l1,l2,l3) 1030 tripV(j,k,i,l2,l3,l1) = tripV(i,j,k,l1,l2,l3) 1031 tripV(k,j,i,l3,l2,l1) = tripV(i,j,k,l1,l2,l3) 1032 tripV(k,i,j,l3,l1,l2) = tripV(i,j,k,l1,l2,l3) 1033C 1034 end do 1035 end do 1036 end do 1037 end do 1038 end do 1039 end do 1040 end if 1041C 1042C Start loop over states 1043C 1044 do index=0, NNM 1045c 1046 if (index.eq.0) then 1047 write(LuOut,9019) 1048 else 1049 write(LuOut,9029) nnm 1050 endif 1051 write(LuOut,*) '' 1052C 1053C array of state indices 1054C 1055 do i=1, NNM 1056 if (index.eq.0) then 1057 istate(i)=0 1058 else if (index.eq.i) then 1059 istate(i)=IEXC 1060 else 1061 istate(i)=0 1062 end if 1063 end do 1064C 1065c INITIALIZE THE WAVEFUNCTIONS 1066C 1067 CALL dfill(nnm*ngrid,zero,wave,1) 1068 CALL dfill(nnm*ngrid,zero,vscf,1) 1069C 1070c CALCULATE THE EIGENVALUES AND EIGENFUNCTIONS OF 1071c DIAGONAL POTENTIAL 1072C 1073 sumE=zero 1074 do mode=1, NNM 1075 do l=1, NGRID 1076 tV(l)=DIAGV(mode,l) 1077 end do 1078C 1079 call collocat(mode,tV,tE,wave,allwave,virtwave,virtE, 1080 * ipvt,rq,dq,xx,a,phi,r,rr,g,v,h,ec,vec,scr,ia, 1081 * gr,temp,twave,nnm,ngrid,ngrid2,nst,istate,index) 1082C 1083 E(mode)=tE 1084 sumE=sumE+tE 1085 end do 1086 ediag(index+1)=sumE 1087 Etot=sumE 1088 emp1=zero 1089 if (ncoup.le.1) goto 50 1090C 1091C START SCF ITERATIONS 1092C 1093 write(LuOut,9050) 1094c 1095 iter=0 1096 1000 continue 1097 iter=iter+1 1098 Eprev=Etot 1099 if (iter.eq.1) Eprev=zero 1100 Etot=zero 1101 sumE=zero 1102 do mode=1, NNM 1103C 1104C effective potential 1105C 1106 call Veffect (mode,vhf,wave,dq,COUPV,TRIPV,NNM,NGRID,ncoup) 1107C 1108 do l=1, NGRID 1109 tV(l)=vhf(l)+diagV(mode,l) 1110 Vscf(mode,l)=vhf(l) 1111 end do 1112C 1113 call collocat(mode,tV,tE,wave,allwave,virtwave,virtE, 1114 * ipvt,rq,dq,xx,a,phi,r,rr,g,v,h,ec,vec,scr,ia,gr, 1115 * temp,twave,nnm,ngrid,ngrid2,nst,istate,index) 1116C 1117 if (tE.lt.zero) then 1118 vcfct1=fct 1119 vcfct=vcfct-1.0D-02 1120 if (vcfct.lt.2.0D-01) then 1121 IF(ga_nodeid().eq.0) 1122 * WRITE(LuOut,*) 'SCALING FACTOR IS LESS THAN 0.2' 1123 CALL errquit('vscfmp: problem with potential',555, 1124 & INPUT_ERR) 1125 end if 1126 if (ga_nodeid().eq.0) write(LuOut,9010) vcfct 1127 goto 500 1128 endif 1129 E(mode)=tE 1130 sumE=sumE + tE 1131 end do 1132C 1133C calculate the SCF correction 1134C 1135 call scfcorr (wave,emp1,dq,coupV,TRIPV,NNM,NGRID,ncoup) 1136C 1137 Etot=sumE-emp1 1138C 1139 write(LuOut,9051) iter,sumE*cm,emp1*cm,Etot*cm 1140c 1141 if (ABS(Etot-Eprev).gt.eps .and. iter.le.mxiter) goto 1000 1142C 1143C exit SCF iterations 1144C 1145 if (ga_nodeid().eq.0) then 1146 if (iter.le.mxiter) then 1147 write(LuOut,*) "VSCF converged in",iter,"iterations" 1148 else 1149 write(LuOut,*) "*** VSCF DID NOT CONVERGE ***" 1150 end if 1151 end if 1152 50 continue 1153c 1154c punch VSCF wavefunctions 1155c 1156 if (ga_nodeid().eq.0.and.print_wv) then 1157 if (index.eq.0) then 1158 write(LuOut,9020) 1159 write(LuOut,9040) vcfct 1160 else 1161 write(LuOut,9030) index 1162 write(LuOut,9040) vcfct 1163 end if 1164 write(LuOut,*) 1165 end if 1166 sumE=ZERO 1167 do mode=1, NNM 1168 if (ga_nodeid().eq.0.and.print_wv) 1169 & write(LuOut,*) mode,istate(mode),E(mode)*cm 1170 do l=1, NGRID 1171 if (ga_nodeid().eq.0.and.print_wv) 1172 & write(LuOut,*) mode,rq(mode,l),wave(mode,l) 1173 end do 1174 sumE=sumE+E(mode) 1175 if (ga_nodeid().eq.0.and.print_wv) write(LuOut,*) 1176 end do 1177c 1178 Etot=sumE-emp1 1179 escf(index+1)=Etot 1180 if (ncoup.le.1) goto 100 1181C 1182C MP2 correction 1183C 1184 call virtstate(index,nnm,nmax,nvst,ivst,ncoup,iexc) 1185C 1186 call mppt(emp2,emp0,virtwave,vmp,virtE,coupV,tripV,Vscf,dq, 1187 * ovrlp,nnm,ngrid,nst,nvst,ivst,istate,jref,jvirt,ncoup) 1188C 1189 Etot=Etot+emp2 1190 emppt(index+1)=Etot 1191 if (ga_nodeid().eq.0) then 1192 write(LuOut,*) 1193 write(LuOut,*) "E(DIAG) :",ediag(index+1)*cm, 1194 & " cm-1 (without mode coupling)" 1195 write(LuOut,*) "E(VSCF) :", (Etot-emp2)*cm, 'cm-1' 1196 if (ncoup.gt.1) write(LuOut,*) "E(VMP2) :", Etot*cm, 'cm-1' 1197 write(LuOut,*) 1198 endif 1199 100 continue 1200C 1201 end do 1202C 1203c calculate vibrational frequencies in cm-1 and print them out 1204c 1205 Nstate=NNM+1 1206 call endiff(freq,frscf,frmp2,ediag,escf,emppt,nnm,Nstate, 1207 * ncoup,iexc) 1208C 1209 RETURN 1210C 1211 9010 FORMAT(/1X,'Scaling coupling potential by ',f4.2) 1212 9019 FORMAT(/1X,'Solving VSCF for the ground vibrational state') 1213 9020 FORMAT(/1X,'Wavefunctions for the ground vibrational state') 1214 9029 FORMAT(/1X,'Solving VSCF for the excited state of mode',I3) 1215 9030 FORMAT(/1X,'Wavefunctions for the excited state of mode',I3) 1216 9040 FORMAT(1X,'(Scaling factor for coupling potential ',f4.2,')') 1217 9050 FORMAT(/1X,'Iteration E(SCF) E(MP1) E(TOTAL)') 1218 9051 FORMAT(4x,i3,5x,f8.2,4x,f8.2,4x,f8.2) 1219 end 1220C 1221c*module vscf *deck collocat 1222 SUBROUTINE collocat(mode,tV,tE,wave,allwave,virtwave,virtE, 1223 * ipvt,x,dx,xx,a,phi,r,rr,g,v,h,e,vec,scr,ia, 1224 * gr,temp,twave,nm,n,n2,nst,istate,index) 1225c 1226 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1227#include "errquit.fh" 1228c 1229#include "global.fh" 1230#include "stdio.fh" 1231C 1232 LOGICAL GOPARR,DSKWRK,MASWRK 1233C 1234 DIMENSION ISTATE(NM) 1235 DIMENSION x(NM,N),dx(NM) 1236 dimension wave(NM,N),allwave(NM+1,NM,N) 1237 dimension tV(N),virtwave(NM,NST,N),virtE(NM,NST) 1238 dimension ipvt(N) 1239 dimension xx(N),A(N),phi(N,N) 1240 dimension R(N,N),RR(N,N),G(N,N),V(N,N) 1241 dimension H(N2),E(N),VEC(N,N) 1242 DIMENSION SCR(N,8),IA(N) 1243 dimension GR(N,N),temp(N),det(2) 1244 dimension twave(N) 1245C 1246C 1247 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, FOUR=4.0D+00) 1248 PARAMETER (TWO=2.0D+00, HALF=0.5D+00, QTR=0.25D+00) 1249 PARAMETER (c=0.7D+00) 1250C 1251c THIS ROUTINE USES THE COLLOCATIONS METHOD TO GENERATE 1252C THE EIGENVALUES AND WAVEFUNCTIONS: E, wave 1253c REF. CPL V153, 1988 pg.98 Yang & Peet 1254c 1255 pi=four*atan(ONE) 1256c 1257c NORMAL COORDINATE ON A GRID 1258c 1259 do l=1, N 1260 xx(l)=x(mode,l) 1261 end do 1262c 1263c GENERATE THE PARAMETERS A'S 1264c 1265 A(1)=(c**2)/((xx(2)-xx(1))**2) 1266 do i=2, N-1 1267 A(i)=(4*(c**2))/((xx(i+1)-xx(i-1))**2) 1268 end do 1269 A(N)=(c**2)/((xx(N)-xx(N-1))**2) 1270c 1271c GENERATE N GAUSSIAN WAVEFUNCTIONS R(i,j) 1272c 1273 do i=1, N 1274 fac=(TWO*A(i)/pi)**QTR 1275 do j=1, N 1276 R(i,j)=fac*exp(-A(i)*(xx(j)-xx(i))**2) 1277 if (abs(R(i,j)).lt.1.0d-99) R(i,j)=zero 1278 RR(i,j)=R(i,j) 1279 end do 1280 end do 1281c 1282c GENERATE THE POTENTIAL MATRIX V(i,j) 1283c 1284 do i=1, N 1285 do j=1, N 1286 V(i,j)=zero 1287 if (j.eq.i) V(i,j)=tV(i) 1288 end do 1289 end do 1290c 1291c GENERATE THE 2ND ORDER DERIVATIVE OF WAVEFUNCTIONS 1292c 1293 do i=1, N 1294 do j=1, N 1295 G(i,j)=zero 1296 dR2=(4*(A(i)**2)*((xx(j)-xx(i))**2))-2*A(i) 1297 G(i,j)=(-HALF)*dR2*R(i,j) 1298 if (abs(G(i,j)).lt.1.0d-99) G(i,j)=zero 1299 end do 1300 end do 1301c 1302c GENERATE THE INVERSE MATRIX OF R(I,J): R-1 1303c 1304 INFO=0 1305 CALL DGEFA(R,N,N,IPVT,INFO) 1306C 1307 IF(INFO.NE.0) THEN 1308 IF(ga_nodeid().eq.0) WRITE(LuOut,*) 'MATRIX R IS SINGULAR' 1309 CALL errquit('collocat: problem generating inverse',555, 1310 & CALC_ERR) 1311 END IF 1312C 1313c options for dgedi routine 1314c job = 11 both determinant and inverse 1315c job = 01 inverse only 1316c job = 10 determinant only 1317c 1318 job=01 1319 call dgedi(R,N,N,ipvt,det,temp,job) 1320c 1321c Multiply kinetic energy matrix G and inverse R 1322c 1323 CALL MRARBR(G,N,N,N,R,N,N,GR,N) 1324 1325c GENERATE THE HAMILTONIAN TO BE DIAGONALIZED : GR-1 + V 1326c (average the off-diagonal terms and put in triangular 1327c form for diagonalization) 1328c 1329 ij=0 1330 do i=1, N 1331 do j=1, i 1332 ij=ij+1 1333 if(i.eq.j) then 1334 H(ij)=GR(i,j)+V(i,j) 1335 VEC(i,j)=H(ij) 1336 else 1337 H(ij)=(GR(i,j)+GR(j,i))/two 1338 VEC(i,j)=H(ij) 1339 VEC(j,i)=H(ij) 1340 endif 1341 end do 1342 end do 1343c 1344c DIAGONALIZE THE PSEDOSPECTRAL MATRIX FOR EIGENVALUES AND EIGENVECTORS 1345c 1346 call util_jacobi(N,VEC,N,E) 1347c 1348c rearrange into descending order 1349c 1350 call eigsort(E,VEC,N,N) 1351c 1352c find wavefunctions 1353c 1354 do i=1, N 1355 do j=1, N 1356 phi(i,j)=VEC(j,i)*RR(j,j) 1357 end do 1358 end do 1359c 1360 do i=1, N 1361 if (i.ge.N-NST+1) then 1362 ist=N-i+1 1363c 1364 area=zero 1365 do j=1, N/2 1366 area=area+phi(i,j) 1367 end do 1368c 1369 if ((mod(ist,2).eq.1).and.(area.lt.zero)) then 1370 do j=1, N 1371 temp(j)=-phi(i,j) 1372 end do 1373 else if ((mod(ist,2).eq.0).and.(area.gt.zero)) then 1374 do j=1, N 1375 temp(j)=-phi(i,j) 1376 end do 1377 else 1378 do j=1, N 1379 temp(j)=phi(i,j) 1380 end do 1381 end if 1382c 1383 call norm (mode,temp,dx,nm,n) 1384c 1385 do j=1, N 1386 virtwave(mode,ist,j)=temp(j) 1387 end do 1388 virtE(mode,ist)=E(i) 1389c 1390 if (i.eq.(N-istate(mode))) then 1391 do j=1, N 1392 twave(j)=temp(j) 1393 end do 1394 tE=E(i) 1395 end if 1396 end if 1397 end do 1398c 1399 do l=1, N 1400 wave(mode,l)=twave(l) 1401 allwave(index+1,mode,l)=twave(l) 1402 end do 1403c 1404 return 1405 end 1406c*module vscf *deck norm 1407 subroutine norm (mode,twave,dx,nm,n) 1408c 1409 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1410C 1411c this routine normalizes the wavefunction obtained from the 1412c collocation grid points for (i) th mode. 1413c 1414 dimension dx(NM) 1415 dimension twave(N) 1416C 1417 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) 1418C 1419C CALCULATE THE NORMALIZATION CONSTANTS USING THE GRID WAVEFUNCTIONS 1420C 1421 wnorm=zero 1422 do m=1, N 1423 temp=dx(mode)*(twave(m)**2) 1424 wnorm=wnorm+temp 1425 end do 1426c 1427c NORMALIZATION OF THE WAVEFUNCTIONS ON THE GRID POINTS 1428c 1429 do m=1, N 1430 twave(m)=(one/sqrt(wnorm))*twave(m) 1431 end do 1432C 1433C CHECK FOR NORMALIZATION: 1434C 1435 area=zero 1436 do l=1, N 1437 area=area+dx(mode)*(twave(l)**2) 1438 end do 1439c 1440 return 1441 end 1442c*module vscf *deck veffect 1443 subroutine Veffect (mode,vscf,wave,dx,Vc,Vtr,nm,n,ncoup) 1444c 1445 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1446c 1447 dimension vscf(N) 1448 dimension wave(NM,N) 1449 dimension dx(NM) 1450 dimension Vc(NM,NM,N,N) 1451 dimension Vtr(NM,NM,NM,N,N,N) 1452c 1453 PARAMETER (ZERO=0.0D+00) 1454c 1455c this routine calculates the scf averaged potential 1456c on the GH-quadrature grid points for (i) th mode. 1457c 1458 do l=1, N 1459 vscf(l)=zero 1460 do j=1, NM 1461 if (j.ne.mode) then 1462 call SCFAVG(sum,wave,dx,Vc,mode,j,l,nm,n) 1463 vscf(l)=vscf(l)+sum 1464 end if 1465 end do 1466 end do 1467 if (ncoup.le.2) goto 100 1468c 1469 do l=1, N 1470 do j=1, NM-1 1471 do k=j+1, NM 1472 if (j.ne.mode .and. k.ne.mode) then 1473 call SCFAVGT(sum,wave,dx,Vtr,mode,j,k,l,nm,n) 1474 vscf(l)=vscf(l)+sum 1475 end if 1476 end do 1477 end do 1478 end do 1479c 1480 100 continue 1481 return 1482 end 1483c*module vscf *deck scfavg 1484 subroutine SCFAVG(sum,wave,dx,Vc,i,j,l,nm,n) 1485c 1486 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1487c 1488c calculate the Vscf(i,l) for mode i and pt. l 1489c 1490 dimension dx(NM) 1491 dimension Vc(NM,NM,N,N) 1492 dimension wave(NM,N) 1493C 1494 PARAMETER (ZERO=0.0D+00) 1495C 1496 sum=zero 1497 do m=1, N 1498 temp1=dx(j)*Vc(i,j,l,m)*(wave(j,m)**2) 1499 sum=sum+temp1 1500 end do 1501C 1502 return 1503 end 1504c*module vscf *deck scfavgt 1505 subroutine SCFAVGT(sum,wave,dx,Vtr,i,j,k,l,nm,n) 1506c 1507 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1508c 1509c calculate the Vscf(i,l) for mode i and pt. l 1510c 1511 dimension dx(NM) 1512 dimension Vtr(NM,NM,NM,N,N,N) 1513 dimension wave(NM,N) 1514C 1515 PARAMETER (ZERO=0.0D+00) 1516C 1517 sum=zero 1518 do m1=1, N 1519 do m2=1, N 1520 temp1=dx(j)*dx(k)*Vtr(i,j,k,l,m1,m2)* 1521 * (wave(j,m1)**2)*(wave(k,m2)**2) 1522 sum=sum+temp1 1523 end do 1524 end do 1525C 1526 return 1527 end 1528c*module vscf *deck scfcorr 1529 subroutine scfcorr (wave,emp1,dx,Vc,Vtr,nm,n,ncoup) 1530c 1531 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1532c 1533c this routine calculates the scf correction 1534c on the GH-quadrature grid points. 1535c 1536 dimension dx(NM) 1537 dimension Vc(NM,NM,N,N) 1538 dimension Vtr(NM,NM,NM,N,N,N) 1539 dimension wave(NM,N) 1540c 1541 PARAMETER (ZERO=0.0D+00) 1542c 1543 emp1=zero 1544 do i=1, NM-1 1545 do j=i+1, NM 1546 do l=1, N 1547 do m=1, N 1548 sum=dx(i)*dx(j)*Vc(i,j,l,m)*(wave(i,l)**2)*(wave(j,m)**2) 1549 emp1=emp1+sum 1550 end do 1551 end do 1552 end do 1553 end do 1554 if (ncoup.le.2) goto 100 1555c 1556 emp1t=zero 1557 do i=1, NM-2 1558 do j=i+1, NM-1 1559 do k=j+1, NM 1560 do l1=1, N 1561 do l2=1, N 1562 do l3=1, N 1563 sum=dx(i)*dx(j)*dx(k)*Vtr(i,j,k,l1,l2,l3)* 1564 * (wave(i,l1)**2)*(wave(j,l2)**2)*(wave(k,l3)**2) 1565 emp1t=emp1t+sum 1566 end do 1567 end do 1568 end do 1569 end do 1570 end do 1571 end do 1572 emp1 = emp1 + emp1t*2.0D+00 1573c 1574 100 continue 1575 return 1576 end 1577c*module vscf *deck virtstate 1578 subroutine virtstate (kstate,nm,nmax,nvst,ivst,ncoup,iexc) 1579c 1580 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1581c 1582 dimension IVST(NVST,NM) 1583c 1584 nost=0 1585c 1586c generate single excitations 1587C 1588 do istate=0, NMAX 1589 do i=1, NM 1590 nost=nost+1 1591 do mode=1, NM 1592 if (mode.eq.i) then 1593 ivst(nost,mode)=istate 1594 else 1595 ivst(nost,mode)=0 1596 end if 1597 end do 1598 if (nost.eq.1) go to 1000 1599 end do 16001000 end do 1601c 1602c generate double excitations 1603C 1604 do jstate=1, NMAX 1605 do istate=1, NMAX 1606c isum=istate+jstate 1607 do i=1, NM-1 1608 do j=i+1, NM 1609 nost=nost+1 1610 do mode=1, NM 1611 if (mode.eq.i) then 1612 ivst(nost,mode)=istate 1613 else if (mode.eq.j) then 1614 ivst(nost,mode)=jstate 1615 else if (mode.eq.kstate) then 1616 ivst(nost,mode)=IEXC 1617 else 1618 ivst(nost,mode)=0 1619 end if 1620 end do 1621 end do 1622 end do 1623 end do 1624 end do 1625 if (ncoup.le.2) goto 2000 1626c 1627c triple excitations if 3body coupling is included 1628C 1629 do istate=1, NMAX 1630 do jstate=1, NMAX 1631 do lstate=1, NMAX 1632c isum=istate+jstate+lstate 1633 do i=1, NM-2 1634 do j=i+1, NM-1 1635 do l=j+1, NM 1636 nost=nost+1 1637 do mode=1, NM 1638 if (mode.eq.i) then 1639 ivst(nost,mode)=istate 1640 else if (mode.eq.j) then 1641 ivst(nost,mode)=jstate 1642 else if (mode.eq.l) then 1643 ivst(nost,mode)=lstate 1644 else if (mode.eq.kstate) then 1645 ivst(nost,mode)=IEXC 1646 else 1647 ivst(nost,mode)=0 1648 end if 1649 end do 1650 end do 1651 end do 1652 end do 1653 end do 1654 end do 1655 end do 1656c 16572000 continue 1658 return 1659 end 1660c*module vscf *deck mppt 1661 subroutine mppt(emp2,emp0,wave,V,E,Vc,Vtr,Vscf,dx,ovrlp, 1662 * nm,n,nst,nvst,ivst,iref,jref,jvirt,ncoup) 1663c 1664c This program calculates the 2nd order MPPT energy correction. 1665c 1666 implicit double precision (a-h,o-z) 1667c 1668#include "global.fh" 1669#include "stdio.fh" 1670c 1671 LOGICAL GOPARR,DSKWRK,MASWRK 1672C 1673 dimension wave(NM,NST,N) 1674 dimension Vc(NM,NM,N,N) 1675 dimension Vtr(NM,NM,NM,N,N,N) 1676 dimension Vscf(NM,N) 1677 dimension E(NM,NST) 1678 dimension dx(NM) 1679 dimension iref(NM) 1680 dimension IVST(NVST,NM) 1681 dimension jvirt(NVST,NM),jref(NM) 1682 dimension emp0(NVST) 1683 dimension v(NVST) 1684 dimension ovrlp(NM) 1685c 1686 PARAMETER (ZERO=0.0D+00, small=1.0D-05, cm=2.19474D+05) 1687c 1688c reference state 1689c 1690 do i=1, NM 1691 iref(i)=iref(i)+1 1692 end do 1693c 1694c compare virtual states with the reference state 1695c 1696 index=0 1697 do j=1, NVST 1698 isum=0 1699 emp0(j)=zero 1700 index=index+1 1701 do i=1, NM 1702 jref(i)=ivst(index,i) 1703 jvirt(j,i)=jref(i)+1 1704 emp0(j)=emp0(j)+E(i,jvirt(j,i)) 1705 if (iref(i).eq.jvirt(j,i)) isum=isum+1 1706 end do 1707 if (isum.eq.NM) istate=j 1708 end do 1709c 1710c generate the overlap between two states i and j 1711c 1712 do j=1, NVST 1713 v(j)=zero 1714 do jmode=1, NM 1715 jref(jmode)=jvirt(j,jmode) 1716 end do 1717 do k=1, NM-1 1718 do l=k+1, NM 1719 call AVG(sum,iref,jref,k,l,dx,wave,Vc,ovrlp,nm,n,nst) 1720 v(j)=v(j)+sum 1721 end do 1722 end do 1723 if (ncoup.gt.2) then 1724 do l1=1, NM-2 1725 do l2=l1+1, NM-1 1726 do l3=l2+1, NM 1727 call AVG3(sum3,iref,jref,l1,l2,l3,dx,wave,Vtr, 1728 * ovrlp,nm,n,nst) 1729 v(j)=v(j)+sum3 1730 end do 1731 end do 1732 end do 1733 end if 1734 do k=1, NM 1735 call AVGHF(sum2,iref,jref,k,dx,wave,Vscf,ovrlp,nm,n,nst) 1736 v(j)=v(j)-sum2 1737 end do 1738 end do 1739c 1740c 1st order energy correction: 1741c 1742 emp1=v(istate) 1743c 1744c 2nd order energy correction: 1745c 1746 nkeep=0 1747 emp2=zero 1748 do j=1, NVST ! state index 1749 if (j.ne.istate) then 1750 sumV=v(j)*v(j) 1751 sumE=emp0(istate)-emp0(j) 1752 if (abs(sumE).lt.small) then 1753 sumV=zero 1754 end if 1755 emp2=emp2+(sumV/sumE) 1756c 1757c to cut down on the computation of coefficients for 1758c 2nd order wavefunction correction, select only the states 1759c which contribute most to the energy correction and use these 1760c states as the basis. 1761c 1762 if (abs(sumV/sumE).gt.1.0D-08) nkeep=nkeep+1 1763 end if 1764 end do 1765 if (ga_nodeid().eq.0) then 1766 write(LuOut,9030) NVST 1767 write(LuOut,9040) nkeep 1768 endif 1769c 1770 return 1771c 1772 9030 format(1x,'Number of virtual states',I6) 1773 9040 format(1x,'Number of selected states',I6) 1774 end 1775c 1776c*module vscf *deck avg 1777 subroutine AVG(sum,iref,jref,k,l,dx,wave,Vc,ovrlp,nm,n,nst) 1778c 1779c this routine calculates the 1780c avg = < Psi | Vc(k,l) | Psi > 1781c 1782 implicit double precision (a-h,o-z) 1783c 1784 dimension iref(NM),jref(NM) 1785 dimension ovrlp(NM) 1786 dimension wave(NM,NST,N) 1787 dimension Vc(NM,NM,N,N) 1788 dimension dx(NM) 1789c 1790 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) 1791c 1792 sum=zero 1793 do kl=1, N 1794 Si=dx(k)*wave(k,iref(k),kl)*wave(k,jref(k),kl) 1795 do ll=1, N 1796 Sj=dx(l)*wave(l,iref(l),ll)*wave(l,jref(l),ll) 1797 sum=sum+Vc(k,l,kl,ll)*Si*Sj 1798 end do 1799 end do 1800c 1801 Sij=one 1802 do m=1, NM 1803 if ((m.eq.k).or.(m.eq.l)) then 1804 ovrlp(m)=one 1805 else 1806 ovrlp(m)=zero 1807 if (iref(m).eq.jref(m)) then 1808 do ll=1, N 1809 ovrlp(m)=ovrlp(m)+ 1810 * dx(m)*wave(m,iref(m),ll)*wave(m,jref(m),ll) 1811 end do 1812 else 1813 Sij=zero 1814 goto 1000 1815 end if 1816 end if 1817 Sij=Sij*ovrlp(m) 1818 end do 18191000 sum=sum*Sij 1820 return 1821 end 1822c*module vscf *deck avg3 1823 subroutine AVG3(sum,iref,jref,l1,l2,l3,dx,wave,Vtr, 1824 * ovrlp,nm,n,nst) 1825c 1826c this routine calculates the 1827c avg = < Psi | Vtrip(l1,l2,l3) | Psi > 1828c 1829 implicit double precision (a-h,o-z) 1830c 1831 dimension iref(NM),jref(NM) 1832 dimension ovrlp(NM) 1833 dimension wave(NM,NST,N) 1834 dimension Vtr(NM,NM,NM,N,N,N) 1835 dimension dx(NM) 1836c 1837 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) 1838c 1839 sum=zero 1840 do ll1=1, N 1841 S1=dx(l1)*wave(l1,iref(l1),ll1)*wave(l1,jref(l1),ll1) 1842 do ll2=1, N 1843 S2=dx(l2)*wave(l2,iref(l2),ll2)*wave(l2,jref(l2),ll2) 1844 do ll3=1, N 1845 S3=dx(l3)*wave(l3,iref(l3),ll3)*wave(l3,jref(l3),ll3) 1846 sum=sum+Vtr(l1,l2,l3,ll1,ll2,ll3)*S1*S2*S3 1847 end do 1848 end do 1849 end do 1850c 1851 Sij=one 1852 do m=1, NM 1853 if ((m.eq.l1).or.(m.eq.l2).or.(m.eq.l3)) then 1854 ovrlp(m)=one 1855 else 1856 ovrlp(m)=zero 1857 if (iref(m).eq.jref(m)) then 1858 do ll=1, N 1859 ovrlp(m)=ovrlp(m)+ 1860 * dx(m)*wave(m,iref(m),ll)*wave(m,jref(m),ll) 1861 end do 1862 else 1863 Sij=zero 1864 goto 1000 1865 end if 1866 end if 1867 Sij=Sij*ovrlp(m) 1868 end do 18691000 sum=sum*Sij 1870 return 1871 end 1872c*module vscf *deck avghf 1873 subroutine AVGHF(sum,iref,jref,k,dx,wave,Vscf,ovrlp,nm,n,nst) 1874c 1875c this routine calculates the 1876c avg = < Psi | Vhf(k) | Psi > 1877c 1878 implicit double precision (a-h,o-z) 1879c 1880 dimension iref(NM),jref(NM) 1881 dimension wave(NM,NST,N) 1882 dimension Vscf(NM,N) 1883 dimension dx(NM) 1884 dimension ovrlp(NM) 1885c 1886 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) 1887c 1888 sum=zero 1889 Sij=one 1890 do m=1, NM 1891 if (m.eq.k) then 1892 ovrlp(m)=one 1893 do ll=1, N 1894 temp=dx(m)*wave(m,jref(m),ll)*wave(m,iref(m),ll) 1895 sum=sum+Vscf(m,ll)*temp 1896 end do 1897 else 1898 ovrlp(m)=zero 1899 if (iref(m).eq.jref(m)) then 1900 do ll=1, N 1901 ovrlp(m)=ovrlp(m)+ 1902 * dx(m)*wave(m,jref(m),ll)*wave(m,iref(m),ll) 1903 end do 1904 else 1905 Sij=zero 1906 goto 1000 1907 end if 1908 end if 1909 Sij=Sij*ovrlp(m) 1910 end do 19111000 sum=sum*Sij 1912 return 1913 end 1914c*module vscf *deck eigsort 1915 SUBROUTINE EIGSORT(D,V,N,NP) 1916c 1917c this routine resorts the eigenvalues from ascending 1918c to descending order, and rearranges the columns of V 1919c correspondingly. 1920c 1921 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1922C 1923 DIMENSION D(NP),V(NP,NP) 1924C 1925 DO 130 I=1,N/2 1926 K=N-I+1 1927 P=D(K) 1928 D(K)=D(I) 1929 D(I)=P 1930 DO 120 J=1,N 1931 P=V(J,K) 1932 V(J,K)=V(J,I) 1933 V(J,I)=P 1934 120 CONTINUE 1935 130 CONTINUE 1936C 1937 RETURN 1938 END 1939c*module vscf *deck endiff 1940 subroutine endiff(freq,frscf,frmp2,diag,emp1,emp2,nm,nstate, 1941 * ncoup,iexc) 1942C 1943 implicit double precision (a-h,o-z) 1944c 1945#include "global.fh" 1946#include "stdio.fh" 1947C 1948 dimension freq(nm), emp1(Nstate), emp2(Nstate), diag(Nstate) 1949 dimension frscf(nm), frmp2(nm) 1950C 1951 PARAMETER (cm=2.19474D+05) 1952C 1953 if (ga_nodeid().eq.0) write(LuOut,*) 1954 if (ga_nodeid().eq.0) 1955 * write(LuOut,*) " RESULTS OF VIBRATIONAL SCF CALCULATION:" 1956 * ," Frequencies, cm-1" 1957 if (ncoup.le.1) then 1958 if (ga_nodeid().eq.0) 1959 * write(LuOut,*) " Mode Harmonic Diagonal" 1960 do imode = 1, NM 1961 hfr = freq(imode)*cm 1962 if (iexc.gt.1) hfr=hfr*iexc 1963 ddiag = (diag(imode+1) - diag(1))*cm 1964 if (ga_nodeid().eq.0) write(LuOut,9010) imode, hfr, ddiag 1965 frscf(imode)=ddiag/cm 1966 frmp2(imode)=frscf(imode) 1967 end do 1968 else 1969 if (ga_nodeid().eq.0) write(LuOut,*) 1970 * " Mode Harmonic Diagonal VSCF", 1971 * " PT2-VSCF" 1972 do imode = 1, NM 1973 hfr = freq(imode)*cm 1974 if (iexc.gt.1) hfr=hfr*iexc 1975 ddiag = (diag(imode+1) - diag(1))*cm 1976 dvscf = (emp1(imode+1) - emp1(1))*cm 1977 demp2 = (emp2(imode+1) - emp2(1))*cm 1978 if (ga_nodeid().eq.0) 1979 * write(LuOut,9020) imode,hfr,ddiag,dvscf,demp2 1980 frscf(imode)=dvscf/cm 1981 frmp2(imode)=demp2/cm 1982 end do 1983 end if 1984C 1985 return 1986C 1987 9010 format (2x,i2,5x,2(f10.2,3x)) 1988 9020 format (2x,i2,5x,4(f10.2,3x)) 1989 end 1990c*module vscf *deck rddiag 1991 SUBROUTINE VSCF_RDDIAG(RQ,DQ,diagV,DMX,DMY,DMZ,NM,N,ndip, 1992 * nstart,ldiag) 1993C 1994 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 1995#include "errquit.fh" 1996c 1997#include "nwc_const.fh" 1998#include "global.fh" 1999#include "mafdecls.fh" 2000#include "stdio.fh" 2001C 2002 DIMENSION DQ(NM),RQ(NM,N),DIAGV(NM,N) 2003 DIMENSION DMX(NM,N),DMY(NM,N),DMZ(NM,N) 2004 dimension nstart(7) 2005C 2006 IF (ga_nodeid().eq.0) THEN 2007 if (ndip.eq.1) then 2008 do i=1, NM 2009 do l=1, N 2010 READ (ldiag,9012,end=9001,err=9001) RQ(i,l), diagV(i,l), 2011 * dmx(i,l),dmy(i,l),dmz(i,l) 2012 nstart(2)=i 2013 nstart(3)=l 2014 enddo 2015 dq(i)=rq(i,2)-rq(i,1) 2016 enddo 2017 else 2018 do i=1, NM 2019 do l=1, N 2020 READ (ldiag,9010,end=9001,err=9001) RQ(i,l), diagV(i,l) 2021 nstart(2)=i 2022 nstart(3)=l 2023 enddo 2024 dq(i)=rq(i,2)-rq(i,1) 2025 enddo 2026 endif 2027 9001 if ((nstart(2).eq.nm).and.(nstart(3).eq.n)) nstart(1)=1 2028 ENDIF 2029c 2030c Now broadcast information to all processors 2031c 2032 master = 0 2033 CALL ga_brdcst(msg_ns1,nstart,7*ma_sizeof(mt_int,1,mt_byte), 2034 * MASTER) 2035 CALL ga_brdcst(msg_rq,RQ,NM*N*ma_sizeof(mt_dbl,1,mt_byte), 2036 * MASTER) 2037 CALL ga_brdcst(msg_dq,DQ,NM*ma_sizeof(mt_dbl,1,mt_byte), 2038 * MASTER) 2039 CALL ga_brdcst(msg_diagv,DIAGV, 2040 * NM*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2041 if (ndip.eq.1) then 2042 CALL ga_brdcst(msg_dmx,DMX, 2043 * NM*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2044 CALL ga_brdcst(msg_dmy,DMY, 2045 * NM*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2046 CALL ga_brdcst(msg_dmz,DMZ, 2047 * NM*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2048 end if 2049C 2050 RETURN 2051 9010 format(5x,f14.8,2x,e14.7) 2052 9012 format(5x,f14.8,2x,e14.7,3(2x,e14.7)) 2053 END 2054c*module vscf *deck rdcoup 2055 SUBROUTINE RDCOUP(coupV,NM,N,NSTART,dm2x,dm2y,dm2z,lcoup,ndip) 2056C 2057 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 2058c 2059#include "global.fh" 2060#include "mafdecls.fh" 2061#include "msgids.fh" 2062C 2063 DIMENSION COUPV(NM,NM,N,N) 2064 DIMENSION dm2x(NM,NM,N,N) 2065 DIMENSION dm2y(NM,NM,N,N) 2066 DIMENSION dm2z(NM,NM,N,N) 2067 dimension nstart(7) 2068C 2069 IF (ga_nodeid().eq.0) THEN 2070 if (ndip.eq.1) then 2071 do i=1, NM 2072 do j=i+1, NM 2073 do l=1, N 2074 do m=1, N 2075 read(lcoup,9022,end=9019,err=9019) coupV(i,j,l,m), 2076 * dm2x(i,j,l,m),dm2y(i,j,l,m),dm2z(i,j,l,m) 2077 nstart(2)=i 2078 nstart(3)=l 2079 nstart(4)=j 2080 nstart(5)=m 2081 enddo 2082 enddo 2083 enddo 2084 enddo 2085 else 2086 do i=1, NM 2087 do j=i+1, NM 2088 do l=1, N 2089 do m=1, N 2090 read(lcoup,9020,end=9019,err=9019) coupV(i,j,l,m) 2091 nstart(2)=i 2092 nstart(3)=l 2093 nstart(4)=j 2094 nstart(5)=m 2095 enddo 2096 enddo 2097 enddo 2098 enddo 2099 end if 2100 9019 if ((nstart(2).eq.nm).and.(nstart(3).eq.n).and.(nstart(4).eq.nm) 2101 & .and.(nstart(5).eq.n)) nstart(1)=2 2102 END IF 2103c 2104c Now broadcast information to all processors 2105c 2106 master = 0 2107 CALL ga_brdcst(msg_ns2,nstart,7*ma_sizeof(mt_int,1,mt_byte), 2108 & MASTER) 2109 CALL ga_brdcst(msg_coupv,COUPV, 2110 * NM*NM*N*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2111 if (ndip.eq.1) then 2112 CALL ga_brdcst(msg_dm2x,dm2x, 2113 * NM*NM*N*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2114 CALL ga_brdcst(msg_dm2y,dm2y, 2115 * NM*NM*N*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2116 CALL ga_brdcst(msg_dm2z,dm2z, 2117 * NM*NM*N*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2118 end if 2119 2120C 2121 RETURN 2122 9020 format(41x,e14.7) 2123 9022 format(41x,e14.7,3(e14.7,2x)) 2124 END 2125c*module vscf *deck rdtrip 2126 SUBROUTINE RDTRIP(tripV,NM,N,NSTART,ltrip) 2127 2128 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 2129 2130#include "global.fh" 2131 2132 integer ltrip 2133 DIMENSION TRIPV(NM,NM,NM,N,N,N) 2134 dimension nstart(7) 2135 2136 IF (ga_nodeid().eq.0) THEN 2137 do i=1, NM 2138 do j=i+1, NM 2139 do k=j+1, NM 2140 do l1=1, N 2141 do l2=1, N 2142 do l3=1, N 2143 read(ltrip,9025,end=9024,err=9024) 2144 & tripV(i,j,k,l1,l2,l3) 2145 nstart(2)=i 2146 nstart(3)=l1 2147 nstart(4)=j 2148 nstart(5)=l2 2149 nstart(6)=k 2150 nstart(7)=l3 2151 enddo 2152 enddo 2153 enddo 2154 enddo 2155 enddo 2156 enddo 2157 9024 if ((nstart(2).eq.nm).and.(nstart(3).eq.n).and.(nstart(4).eq.nm) 2158 & .and.(nstart(5).eq.n).and.(nstart(6).eq.nm) 2159 & .and.(nstart(7).eq.n)) nstart(1)=3 2160 END IF 2161c 2162c Now broadcast information to all processors 2163c 2164 master = 0 2165 CALL ga_brdcst(msg_ns3,nstart,7*ma_sizeof(mt_int,1,mt_byte), 2166 & MASTER) 2167 CALL ga_brdcst(msg_tripv,TRIPV, 2168 * NM*NM*NM*N*N*N*ma_sizeof(mt_dbl,1,mt_byte),MASTER) 2169 RETURN 2170 9025 format(61x,e14.7) 2171 END 2172c*module vscf *deck intens 2173 SUBROUTINE INTENS(SINTSCF,UX,UY,UZ,FRSCF,FRMP2,DQ,WAVE,NM,N) 2174C 2175 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 2176c 2177#include "global.fh" 2178#include "stdio.fh" 2179C 2180c THIS ROUTINE CALCULATES IR intensity for i-th MODE 2181c USING THE following equation: 2182c 2183c Intensity(i) = |<ground psi(i)| u(i) |excited psi(i)>|^2 2184c 2185c where psi(i) = wavefunction corresponding to mode i 2186c u(i) = dipole moment along normal coordinate i 2187c 2188 dimension dQ(NM) 2189 dimension wave(NM+1,NM,N) ! VSCF wavefunctions 2190 dimension frscf(NM) ! VSCF frequencies 2191 dimension frmp2(NM) ! VMP2 frequencies 2192 dimension ux(NM,N),uy(NM,N),uz(NM,N) ! dipole moment components 2193 dimension sIntscf(NM) ! VSCF Intensity 2194c 2195 PARAMETER (ZERO=0.0D+00) 2196 PARAMETER (const=2.5048D+00, cm=2.19474D+05) 2197c 2198c const = 8*pi^3*Navogadro*(E-18)^2*E-5/(3*h*c) in CGS 2199c 2200 if (ga_nodeid().eq.0) then 2201 write(LuOut,*) 2202 write(LuOut,*) " IR INTENSITIES CALCULATED USING DIPOLE MOMENT:" 2203 write(LuOut,*) " mode frequency, cm-1 intensity, km/mol" 2204 endif 2205c 2206c Note: dipole assumed to be in debeye 2207c 2208 do j=1, NM 2209 sIntx = zero 2210 sInty = zero 2211 sIntz = zero 2212 do l=1, N 2213 sIntx = sIntx + 2214 * ux(j,l)*wave(1,j,l)*wave(j+1,j,l)*dQ(j) 2215 sInty = sInty + 2216 * uy(j,l)*wave(1,j,l)*wave(j+1,j,l)*dQ(j) 2217 sIntz = sIntz + 2218 * uz(j,l)*wave(1,j,l)*wave(j+1,j,l)*dQ(j) 2219 end do 2220 sIntscf(j) = sIntx*sIntx+sInty*sInty+sIntz*sIntz 2221 sIntscf(j) = frscf(j)*cm*const*sIntscf(j) 2222 if (ga_nodeid().eq.0) 2223 * write(LuOut,9000) j, frmp2(j)*cm, sIntscf(j) 2224 end do 2225c 2226 if (ga_nodeid().eq.0) then 2227 write(LuOut,*) 2228 write(LuOut,*) 2229 * " (Note: SCF dipole moments used to calculate intensities)" 2230 endif 2231 return 2232c 2233 9000 format (2x,i2,5x,f10.2,4x,f10.2) 2234 end 2235c*module vscf *deck dintens 2236 SUBROUTINE DINTENS(SINTSCF,DDM,DDER,FRSCF,FRMP2, 2237 * VEC,Q,DQ,WAVE,NM,N,NC) 2238C 2239 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 2240c 2241#include "global.fh" 2242#include "stdio.fh" 2243C 2244c THIS ROUTINE CALCULATES IR intensity for i-th MODE 2245c USING THE following equation: 2246c 2247c Intensity(i) = ( du/dQ(i) |<gr. psi(i)| Q(i) |ex. psi(i)>| )^2 2248c 2249c where du/dQ(i) = dipole derivative at equilibrium 2250c 2251 dimension dQ(NM) 2252 dimension Q(NM,N) 2253 dimension VEC(NC,NC) 2254 dimension wave(NM+1,NM,N) ! VSCF wavefunctions 2255 dimension frscf(NM) ! VSCF frequencies 2256 dimension frmp2(NM) ! VMP2 frequencies 2257 dimension sIntscf(NM) ! VSCF Intensity 2258 dimension ddm(NC*3) ! dipole derivative matrix in cartezians 2259 dimension dder(NM) ! square of dipole moment derivative 2260c 2261 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00) 2262 PARAMETER (AMU=1.8229D+03, cm=2.19474D+05) 2263 PARAMETER (const=2.5048D+00, bohr=5.2918D-01) 2264c 2265c const = 8*pi^3*Navogadro*(E-18)^2*E-5/(3*h*c) in CGS 2266c 2267 FACT = one/AMU 2268 NSTART=NC-NM+1 2269C 2270 if (ga_nodeid().eq.0) write(LuOut,*) 2271 if (ga_nodeid().eq.0) write(LuOut,*) 2272 * " IR INTENSITIES CALCULATED USING DIPOLE DERIVATIVE:" 2273c 2274 do i = NSTART,NC 2275 im=nc-i+1 2276 DDX = DDOT(NC,VEC(1,i),1,DDM(1),3) 2277 DDY = DDOT(NC,VEC(1,i),1,DDM(2),3) 2278 DDZ = DDOT(NC,VEC(1,i),1,DDM(3),3) 2279 dder(im) = DDX*DDX + DDY*DDY + DDZ*DDZ 2280 enddo 2281c 2282 if (ga_nodeid().eq.0) write(LuOut,*) 2283 * " mode frequency, cm-1 intensity, km/mol" 2284c 2285 do j=1, NM 2286 sIntscf(j) = zero 2287 do l=1, N 2288 sIntscf(j) = sIntscf(j) + 2289 * Q(j,l)*wave(1,j,l)*wave(j+1,j,l)*dQ(j) 2290 end do 2291c 2292 sIntscf(j) = frscf(j)*cm*sIntscf(j)**2 2293 sIntscf(j) = sIntscf(j)*dder(j)*const*fact*bohr*bohr 2294 if (ga_nodeid().eq.0) 2295 * write(LuOut,9000) j, frmp2(j)*cm, sIntscf(j) 2296c 2297 end do 2298c 2299 return 2300c 2301 9000 format (2x,i2,5x,f10.2,4x,f10.2) 2302 end 2303C*MODULE MTHLIB *DECK DGEFA 2304 SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) 2305 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 2306 DIMENSION A(LDA,*),IPVT(*) 2307C 2308C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. 2309C 2310C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED 2311C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. 2312C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) . 2313C 2314C ON ENTRY 2315C 2316C A DOUBLE PRECISION(LDA, N) 2317C THE MATRIX TO BE FACTORED. 2318C 2319C LDA INTEGER 2320C THE LEADING DIMENSION OF THE ARRAY A . 2321C 2322C N INTEGER 2323C THE ORDER OF THE MATRIX A . 2324C 2325C ON RETURN 2326C 2327C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS 2328C WHICH WERE USED TO OBTAIN IT. 2329C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE 2330C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER 2331C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. 2332C 2333C IPVT INTEGER(N) 2334C AN INTEGER VECTOR OF PIVOT INDICES. 2335C 2336C INFO INTEGER 2337C = 0 NORMAL VALUE. 2338C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR 2339C CONDITION FOR THIS SUBROUTINE, BUT IT DOES 2340C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO 2341C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE 2342C INDICATION OF SINGULARITY. 2343C 2344C LINPACK. THIS VERSION DATED 08/14/78 . 2345C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. 2346C 2347C SUBROUTINES AND FUNCTIONS 2348C 2349C BLAS DAXPY,DSCAL,IDAMAX 2350C 2351C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 2352C 2353 INFO = 0 2354 NM1 = N - 1 2355 IF (NM1 .LT. 1) GO TO 70 2356 DO 60 K = 1, NM1 2357 KP1 = K + 1 2358C 2359C FIND L = PIVOT INDEX 2360C 2361 L = IDAMAX(N-K+1,A(K,K),1) + K - 1 2362 IPVT(K) = L 2363C 2364C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED 2365C 2366 IF (A(L,K) .EQ. 0.0D+00) GO TO 40 2367C 2368C INTERCHANGE IF NECESSARY 2369C 2370 IF (L .EQ. K) GO TO 10 2371 T = A(L,K) 2372 A(L,K) = A(K,K) 2373 A(K,K) = T 2374 10 CONTINUE 2375C 2376C COMPUTE MULTIPLIERS 2377C 2378 T = -1.0D+00/A(K,K) 2379 CALL DSCAL(N-K,T,A(K+1,K),1) 2380C 2381C ROW ELIMINATION WITH COLUMN INDEXING 2382C 2383 DO 30 J = KP1, N 2384 T = A(L,J) 2385 IF (L .EQ. K) GO TO 20 2386 A(L,J) = A(K,J) 2387 A(K,J) = T 2388 20 CONTINUE 2389 CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 2390 30 CONTINUE 2391 GO TO 50 2392 40 CONTINUE 2393 INFO = K 2394 50 CONTINUE 2395 60 CONTINUE 2396 70 CONTINUE 2397 IPVT(N) = N 2398 IF (A(N,N) .EQ. 0.0D+00) INFO = N 2399 RETURN 2400 END 2401C*MODULE MTHLIB *DECK DGEDI 2402 SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB) 2403 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 2404 DIMENSION A(LDA,*),DET(2),WORK(*),IPVT(*) 2405C 2406C DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX 2407C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. 2408C 2409C ON ENTRY 2410C 2411C A DOUBLE PRECISION(LDA, N) 2412C THE OUTPUT FROM DGECO OR DGEFA. 2413C 2414C LDA INTEGER 2415C THE LEADING DIMENSION OF THE ARRAY A . 2416C 2417C N INTEGER 2418C THE ORDER OF THE MATRIX A . 2419C 2420C IPVT INTEGER(N) 2421C THE PIVOT VECTOR FROM DGECO OR DGEFA. 2422C 2423C WORK DOUBLE PRECISION(N) 2424C WORK VECTOR. CONTENTS DESTROYED. 2425C 2426C JOB INTEGER 2427C = 11 BOTH DETERMINANT AND INVERSE. 2428C = 01 INVERSE ONLY. 2429C = 10 DETERMINANT ONLY. 2430C 2431C ON RETURN 2432C 2433C A INVERSE OF ORIGINAL MATRIX IF REQUESTED. 2434C OTHERWISE UNCHANGED. 2435C 2436C DET DOUBLE PRECISION(2) 2437C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. 2438C OTHERWISE NOT REFERENCED. 2439C DETERMINANT = DET(1) * 10.0**DET(2) 2440C WITH 1.0 .LE. ABS(DET(1)) .LT. 10.0 2441C OR DET(1) .EQ. 0.0 . 2442C 2443C ERROR CONDITION 2444C 2445C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS 2446C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. 2447C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY 2448C AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET 2449C INFO .EQ. 0 . 2450C 2451C LINPACK. THIS VERSION DATED 08/14/78 . 2452C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. 2453C 2454C SUBROUTINES AND FUNCTIONS 2455C 2456C BLAS DAXPY,DSCAL,DSWAP 2457C FORTRAN ABS,MOD 2458C 2459C COMPUTE DETERMINANT 2460C 2461 IF (JOB/10 .EQ. 0) GO TO 70 2462 DET(1) = 1.0D+00 2463 DET(2) = 0.0D+00 2464 TEN = 10.0D+00 2465 DO 50 I = 1, N 2466 IF (IPVT(I) .NE. I) DET(1) = -DET(1) 2467 DET(1) = A(I,I)*DET(1) 2468C ...EXIT 2469 IF (DET(1) .EQ. 0.0D+00) GO TO 60 2470 10 IF (ABS(DET(1)) .GE. 1.0D+00) GO TO 20 2471 DET(1) = TEN*DET(1) 2472 DET(2) = DET(2) - 1.0D+00 2473 GO TO 10 2474 20 CONTINUE 2475 30 IF (ABS(DET(1)) .LT. TEN) GO TO 40 2476 DET(1) = DET(1)/TEN 2477 DET(2) = DET(2) + 1.0D+00 2478 GO TO 30 2479 40 CONTINUE 2480 50 CONTINUE 2481 60 CONTINUE 2482 70 CONTINUE 2483C 2484C COMPUTE INVERSE(U) 2485C 2486 IF (MOD(JOB,10) .EQ. 0) GO TO 150 2487 DO 100 K = 1, N 2488 A(K,K) = 1.0D+00/A(K,K) 2489 T = -A(K,K) 2490 CALL DSCAL(K-1,T,A(1,K),1) 2491 KP1 = K + 1 2492 IF (N .LT. KP1) GO TO 90 2493 DO 80 J = KP1, N 2494 T = A(K,J) 2495 A(K,J) = 0.0D+00 2496 CALL DAXPY(K,T,A(1,K),1,A(1,J),1) 2497 80 CONTINUE 2498 90 CONTINUE 2499 100 CONTINUE 2500C 2501C FORM INVERSE(U)*INVERSE(L) 2502C 2503 NM1 = N - 1 2504 IF (NM1 .LT. 1) GO TO 140 2505 DO 130 KB = 1, NM1 2506 K = N - KB 2507 KP1 = K + 1 2508 DO 110 I = KP1, N 2509 WORK(I) = A(I,K) 2510 A(I,K) = 0.0D+00 2511 110 CONTINUE 2512 DO 120 J = KP1, N 2513 T = WORK(J) 2514 CALL DAXPY(N,T,A(1,J),1,A(1,K),1) 2515 120 CONTINUE 2516 L = IPVT(K) 2517 IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1) 2518 130 CONTINUE 2519 140 CONTINUE 2520 150 CONTINUE 2521 RETURN 2522 END 2523C*MODULE MTHLIB *DECK MRARBR 2524 SUBROUTINE MRARBR(A,LDA,NA,MA,B,LDB,MB,C,LDC) 2525C 2526 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 2527#include "errquit.fh" 2528C 2529#include "global.fh" 2530c 2531 DIMENSION A(LDA,MA),B(LDB,MB),C(LDC,MB) 2532C 2533 PARAMETER (ZERO=0.0D+00) 2534*BL3 PARAMETER (ONE=1.0D+00) 2535C 2536C ----- SEQUENTIAL MATRIX MULTIPLY: C = A * B ----- 2537C ----- SEE ALSO -PRARBR- ROUTINE ----- 2538C A - THE INPUT REAL NA X MA RECTANGULAR MATRIX 2539C B - THE INPUT REAL MA X MB RECTANGULAR MATRIX 2540C C - THE OUTPUT PRODUCT NA X MB MATRIX 2541C LDA - THE LEADING DIMENSION OF MATRIX A 2542C NA - THE EFFECTIVE ROW DIMENSION OF MATRICES A AND C 2543C MA - THE COLUMN DIMENSION OF MATRIX A, 2544C AND EFFECTIVE ROW DIMENSION OF MATRIX B 2545C LDB - THE LEADING DIMENSION OF MATRIX B 2546C MB - THE COLUMN DIMENSION OF MATRICES B AND C 2547C LDC - THE LEADING DIMENSION OF MATRIX C 2548C AUTHOR = STEVE ELBERT, 31 OCT 1979 2549C 2550 IF(LDA.LT.NA .OR. LDB.LT.MA .OR. LDC.LT.NA) GO TO 800 2551C 2552*BL3 CALL DGEMM('N','N',NA,MB,MA,ONE,A,LDA,B,LDB,ZERO,C,LDC) 2553*BL3 IF(ONE.NE.ZERO) RETURN 2554C 2555 M=MA 2556 IF(MOD(M,2).NE.0) M=M-1 2557C 2558 DO 300 I=1,NA 2559 DO 200 J=1,MB 2560 CIJ=ZERO 2561 IF(M.NE.MA) CIJ=A(I,MA)*B(MA,J) 2562 IF(MA.GT.1) THEN 2563 DO 100 K=1,M,2 2564 CIJ=CIJ + A(I,K)*B(K,J) + A(I,K+1)*B(K+1,J) 2565 100 CONTINUE 2566 END IF 2567 C(I,J)=CIJ 2568 200 CONTINUE 2569 300 CONTINUE 2570 RETURN 2571C 2572 800 CONTINUE 2573 IF(ga_nodeid().eq.0) WRITE(6,900) LDA,NA,MA,LDB,MB,LDC 2574 CALL errquit("mrarbr: bad value of a leading dimension",555, 2575 & INPUT_ERR) 2576 STOP 2577C 2578 900 FORMAT(/1X,'ERROR IN CALL TO -MRARBR-'/ 2579 * 1X,'LDA,NA,MA,LDB,MB,LDC=',6I10) 2580 END 2581c 2582 subroutine vscf_restart(rtdb,restart,ncoup,nstart,ldiag,lcoup, 2583 & ltrip,nnm,ngrid,rq,dq,diagv,dm1x,dm1y, 2584 & dm1z,dm2x,dm2y,dm2z,coupv,tripv,ndip) 2585 implicit double precision(a-h,o-z) 2586#include "errquit.fh" 2587#include "nwc_const.fh" 2588#include "geomP.fh" 2589#include "geom.fh" 2590#include "global.fh" 2591#include "mafdecls.fh" 2592#include "stdio.fh" 2593#include "rtdb.fh" 2594c 2595 dimension rq(nnm,ngrid),dq(nnm),diagv(nnm,ngrid) 2596 dimension dm1x(nnm,ngrid),dm1y(nnm,ngrid),dm1z(nnm,ngrid) 2597 dimension dm2x(nnm,nnm,ngrid,ngrid),dm2y(nnm,nnm,ngrid,ngrid) 2598 dimension dm2z(nnm,nnm,ngrid,ngrid),coupv(nnm,nnm,ngrid,ngrid) 2599 dimension tripv(nnm,nnm,nnm,ngrid,ngrid,ngrid) 2600c 2601 character*255 namediag, namecoup, nametrip 2602 dimension nstart(7) 2603 logical restart,diag_exists,coup_exists,trip_exists 2604 integer rtdb 2605c 2606 do ii=1,7 2607 nstart(ii) = -1 2608 enddo 2609c 2610c Open files for diagonal, coupling, and triple potentials 2611c 2612 ldiag = 81 2613 lcoup = 82 2614 ltrip = 83 2615c 2616 if (ga_nodeid().eq.0) then 2617 diag_exists=.false. 2618 coup_exists=.false. 2619 trip_exists=.false. 2620 if (ncoup.ge.1) then 2621 call util_file_name('diag',.false.,.false.,namediag) 2622 inquire(file=namediag,exist=diag_exists) 2623 OPEN (UNIT=ldiag, FORM='FORMATTED', FILE=namediag, 2624 * ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 2625 REWIND (UNIT=ldiag) 2626 endif 2627 if (ncoup.ge.2) then 2628 call util_file_name('coup',.false.,.false.,namecoup) 2629 inquire(file=namecoup,exist=coup_exists) 2630 OPEN (UNIT=lcoup, FORM='FORMATTED', FILE=namecoup, 2631 * ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 2632 REWIND (UNIT=lcoup) 2633 endif 2634 if (ncoup.ge.3) then 2635 call util_file_name('trip',.false.,.false.,nametrip) 2636 inquire(file=nametrip,exist=trip_exists) 2637 OPEN (UNIT=ltrip, FORM='FORMATTED', FILE=nametrip, 2638 * ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 2639 REWIND (UNIT=ltrip) 2640 endif 2641 end if 2642c 2643c If there is restart data, make sure all the processors 2644c are aware of it and get it from the files. 2645c 2646 master=0 2647 CALL ga_brdcst(msg_diagv,diag_exists, 2648 * ma_sizeof(mt_log,1,mt_byte),MASTER) 2649 CALL ga_brdcst(msg_diagv,coup_exists, 2650 * ma_sizeof(mt_log,1,mt_byte),MASTER) 2651 CALL ga_brdcst(msg_diagv,trip_exists, 2652 * ma_sizeof(mt_log,1,mt_byte),MASTER) 2653 if (diag_exists) then 2654 if(ga_nodeid().eq.0) write(LuOut,9030) 2655 CALL VSCF_RDDIAG(rq,dq,diagv,dm1x,dm1y,dm1z,nnm,ngrid,ndip, 2656 * nstart,ldiag) 2657 endif 2658 if (coup_exists.and.ncoup.ge.2) then 2659 if(ga_nodeid().eq.0) write(LuOut,9040) 2660 CALL RDCOUP(coupv,nnm,ngrid,nstart,dm2x,dm2y,dm2z,lcoup,ndip) 2661 endif 2662 if (trip_exists.and.ncoup.eq.3) then 2663 if(ga_nodeid().eq.0) write(LuOut,9050) 2664 CALL RDTRIP(TRIPV,NNM,NGRID,nstart,ltrip) 2665 endif 2666c 2667c Restart if user specified restart keyword and we are reusing 2668c the old rtdb 2669c 2670 if (.not. rtdb_get(rtdb,'vscf:restart',mt_log,1,restart)) 2671 & restart=.false. 2672c 2673 return 2674 9030 format(/1x,'Reading diagonal energy data for restart...') 2675 9040 format(/1x,'Reading coupling energy data for restart...') 2676 9050 format(/1x,'Reading 3-body energy data for restart...') 2677 end 2678