1 subroutine ccsd_imaginary(d_a0,d_f1,d_v2,d_d1, 2 1 d_t1,d_t2,d_lambda1,d_lambda2,d_tr1,d_tr2, 3 2 k_a0_offset,k_f1_offset,k_v2_offset,k_d1_offset, 4 4 k_t1_offset,k_t2_offset,k_l1_offset,k_l2_offset, 5 6 k_tr1_offset,k_tr2_offset, 6 7 size_tr1,size_tr2) 7 implicit none 8#include "global.fh" 9#include "mafdecls.fh" 10#include "util.fh" 11#include "errquit.fh" 12#include "stdio.fh" 13#include "tce.fh" 14#include "tce_main.fh" 15#include "tce_prop.fh" 16#include "tce_restart.fh" 17c 18 integer i,j,dummy,axis 19 integer omegacount 20 integer omegasign 21 integer dynfreq 22 integer dynaxis 23 integer irrep_g 24 parameter (irrep_g=0) 25 integer d_a0,d_f1,d_v2,d_d1(3) 26 integer d_t1,d_t2,d_lambda1,d_lambda2 27 integer d_tr1(9),d_tr2(9) 28 integer k_a0_offset,k_f1_offset,k_v2_offset,k_d1_offset(3) 29 integer k_t1_offset,k_t2_offset,k_l1_offset,k_l2_offset 30 integer k_tr1_offset(3),k_tr2_offset(3) 31 integer size_tr1(3),size_tr2(3) 32 integer sym_abelian_axis 33 integer zero,one 34 external sym_abelian_axis 35 double precision pi 36 parameter(pi = 3.14159265358979323846264338327950288419D0) 37 parameter(zero = 0) 38 parameter(one = 0) 39 double precision omega 40 double precision threeoverpi 41 parameter (threeoverpi=3.0d0/pi) 42 logical nodezero,guess 43 character*4 irrepname 44 character*3 axisname(3) ! Axis 45 data axisname/'X','Y','Z'/ 46c 47 nodezero=(ga_nodeid().eq.0) 48c 49 if(nodezero) then 50 write(6,*) 51 write(6,*) 'Casimir-Polder integration points are:' 52 do i = 1, inumfreq 53 write(6,999) i,ifreq(i) 54 enddo 55 endif 56c 57 do omegacount=1,inumfreq 58 omega = ifreq(omegacount) 59 do axis = 1, 3 60 if (respaxis(axis)) then 61 irrep_d=sym_abelian_axis(geom,axis) 62 call sym_irrepname(geom,irrep_d+1,irrepname) 63 if (nodezero.and.util_print('mod1',print_default)) then 64 write(LuOut,*) 65 write(LuOut,9440) axisname(axis),irrepname 66 endif 67 irrep_o=irrep_d 68 irrep_x=irrep_d 69c 70 if (nodezero) write(LuOut,9431) omega 71c 72 if (cc_ir_alg.eq.1) then 73c 74c REAL COMPONENT 75c 76 dynaxis = 6 77 if (guess_ir_real.and.(omegacount.eq.1)) then 78 call ccsd_ir_guess(d_d1,d_t1,d_t2,d_tr1,d_tr2, 79 1 k_d1_offset,k_t1_offset,k_t2_offset, 80 2 k_tr1_offset,k_tr2_offset, 81 4 size_tr1,size_tr2,zero,axis,omega) 82 endif ! guess_ir_real 83 call ccsd_ir_real_sq_it(axis,dynaxis,omega, 84 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 85 2 k_f1_offset,k_v2_offset,k_d1_offset, 86 3 k_t1_offset,k_t2_offset, 87 4 k_tr1_offset,k_tr2_offset, 88 5 size_tr1,size_tr2) 89c 90c IMAGINARY COMPONENT 91c 92 if (guess_ir_imag.and.(omegacount.eq.1)) then 93 call ccsd_ir_get_imag_from_real(axis,omega, 94 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 95 2 k_f1_offset,k_v2_offset,k_d1_offset, 96 3 k_t1_offset,k_t2_offset,k_tr1_offset,k_tr2_offset, 97 4 size_tr1,size_tr2) 98 endif 99c 100 dynaxis = 0 ! positive omega_I 101 call ccsd_ir_imag_it(axis,dynaxis,omega, 102 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 103 2 k_f1_offset,k_v2_offset,k_d1_offset, 104 3 k_t1_offset,k_t2_offset, 105 4 k_tr1_offset,k_tr2_offset, 106 5 size_tr1,size_tr2) 107 call daxfile(1,-1.0d0,d_tr1(axis+0),d_tr1(axis+3), 108 1 size_tr1(axis)) 109 call daxfile(1,-1.0d0,d_tr2(axis+0),d_tr2(axis+3), 110 2 size_tr2(axis)) 111#ifdef DEBUG_ONLY 112 dynaxis = 3 ! is for negative omega_I 113 call ccsd_ir_imag_it(axis,dynaxis,(-1.0d0*omega), 114 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 115 2 k_f1_offset,k_v2_offset,k_d1_offset, 116 3 k_t1_offset,k_t2_offset, 117 4 k_tr1_offset,k_tr2_offset, 118 5 size_tr1,size_tr2) 119#endif 120c 121 elseif (cc_ir_alg.eq.2) then 122 if (omega.gt.(1.0d-3)) then 123c 124c IMAGINARY COMPONENT 125c 126 dynaxis = 0 ! positive omega_I 127 call ccsd_ir_imag_sq_it(axis,dynaxis,1.0d0/omega, 128 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 129 2 k_f1_offset,k_v2_offset,k_d1_offset, 130 3 k_t1_offset,k_t2_offset, 131 4 k_tr1_offset,k_tr2_offset, 132 5 size_tr1,size_tr2) 133 call daxfile(1,-1.0d0,d_tr1(axis+0),d_tr1(axis+3), 134 1 size_tr1(axis)) 135 call daxfile(1,-1.0d0,d_tr2(axis+0),d_tr2(axis+3), 136 2 size_tr2(axis)) 137#ifdef DEBUG_ONLY 138 dynaxis = 3 ! negative omega_I 139 call ccsd_ir_imag_sq_it(axis,dynaxis,(-1.0d0/omega), 140 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 141 2 k_f1_offset,k_v2_offset,k_d1_offset, 142 3 k_t1_offset,k_t2_offset, 143 4 k_tr1_offset,k_tr2_offset, 144 5 size_tr1,size_tr2) 145#endif 146c 147c REAL COMPONENT 148c 149 dynaxis = 6 ! offset to the third set of response for real component 150 call ccsd_ir_real_it(axis,dynaxis,1.0d0/omega, 151 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 152 2 k_f1_offset,k_v2_offset,k_d1_offset, 153 3 k_t1_offset,k_t2_offset, 154 4 k_tr1_offset,k_tr2_offset, 155 5 size_tr1,size_tr2) 156 elseif (omega.eq.(0.0d0)) then 157c 158c IMAGINARY COMPONENT 159c 160 call tce_zero(d_tr1(axis+0),size_tr1(axis)) 161 call tce_zero(d_tr2(axis+0),size_tr2(axis)) 162 call tce_zero(d_tr1(axis+3),size_tr1(axis)) 163 call tce_zero(d_tr2(axis+3),size_tr2(axis)) 164c 165c REAL COMPONENT 166c 167 dynaxis = 6 ! offset to the third set of response for real component 168 call ccsd_lr_iter(axis,dynaxis,omega, 169 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 170 2 k_f1_offset,k_v2_offset,k_d1_offset, 171 3 k_t1_offset,k_t2_offset,k_tr1_offset,k_tr2_offset, 172 4 size_tr1,size_tr2) 173 else 174 if (nodezero) then 175 write(9121) 'Frequency is too close to zero.' 176 endif 177 endif 178c 179 elseif (cc_ir_alg.eq.3) then 180c 181 guess = .true. !(omegacount.le.2) 182c 183 call ccsd_ir_coupled_it(guess,axis,omega, 184 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 185 2 k_f1_offset,k_v2_offset,k_d1_offset, 186 3 k_t1_offset,k_t2_offset, 187 4 k_tr1_offset,k_tr2_offset, 188 5 size_tr1,size_tr2) 189c 190 endif ! cc_ir_alg 191c 192 endif ! respaxis(axis) 193 enddo ! axis loop 194c 195c CCSD-IR evaluation step 196c 197 call ccsd_ir_eval(omegacount,d_a0,d_f1,d_v2,d_d1, 198 1 d_t1,d_t2,d_lambda1,d_lambda2, 199 2 d_tr1,d_tr2,k_a0_offset, 200 3 k_f1_offset,k_v2_offset,k_d1_offset, 201 4 k_t1_offset,k_t2_offset, 202 5 k_l1_offset,k_l2_offset, 203 6 k_tr1_offset,k_tr2_offset) 204c 205 enddo ! omegacount loop 206c 207c Casimir-Polder integration for determination of C6 coefficient 208c 209 if (ifreqauto) then 210 do i = 1, 8 211 integral(i) = 0.0d0 212 do omegacount = 1,inumfreq-1 213 num0 = ifreqval(omegacount,i) 214 num1 = ifreqval(omegacount+1,i) 215 den0 = ifreq(omegacount) 216 den1 = ifreq(omegacount+1) 217 num0sq = num0*num0 218 num1sq = num1*num1 219 numtot = num1sq+num0sq 220 dentot = den1-den0 221 integral(i) = integral(i) + 0.5d0*numtot*dentot 222 enddo 223 integral(i) = integral(i)*threeoverpi 224 enddo 225c 226 if (nodezero) write(LuOut,9435) "CCSD Imaginary Response", 227 1 integral(1),integral(2),integral(3),integral(4), 228 2 integral(5),integral(6),integral(7),integral(8) 229c 230 endif 231c 232 999 format(3x,'ifreq(',i2,') = ',f14.8) 233 9020 format(1x,'Cpu & wall time / sec',2f15.1) 234 9100 format(1x,i4,2f18.13,2f8.1) 235 9120 format(1x,A) 236 9121 format(/,1x,A) 237 9122 format(1x,A,i4) 238 9420 format(1x,i4,f25.13,2f8.1) 239 9431 format(/,1x,'Frequency = ',f15.7,' / au') 240 9435 format(/,1x,A,' C6 coefficients ',/ 241 1 1x,'--------------------------------',/ 242 2 1x,'C6(XX) ',f15.7,/ 243 3 1x,'C6(YY) ',f15.7,/ 244 4 1x,'C6(ZZ) ',f15.7,/ 245 5 1x,'C6(XY) ',f15.7,/ 246 6 1x,'C6(XZ) ',f15.7,/ 247 7 1x,'C6(YZ) ',f15.7,/ 248 8 1x,'C6(AVG) ',f15.7,/ 249 9 1x,'C6(ANI) ',f15.7,/ 250 1 1x,'--------------------------------') 251 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 252 return 253 end 254 255 256 subroutine ccsd_ir_real_sq_it(axis,dynaxis,omega, 257 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 258 2 k_f1_offset,k_v2_offset,k_d1_offset, 259 3 k_t1_offset,k_t2_offset,k_tr1_offset,k_tr2_offset, 260 4 size_tr1,size_tr2) 261 implicit none 262#include "global.fh" 263#include "mafdecls.fh" 264#include "util.fh" 265#include "errquit.fh" 266#include "stdio.fh" 267#include "tce.fh" 268#include "tce_main.fh" 269#include "tce_diis.fh" 270#include "tce_prop.fh" 271#include "tce_restart.fh" 272c 273 integer i,j,dummy,axis 274 integer omegacount 275 integer omegasign 276 integer dynfreq 277 integer dynaxis 278 integer axisA 279 integer axisB 280 integer irrep_g 281 parameter (irrep_g=0) 282 integer d_f1,d_v2,d_d1(3),d_t1,d_t2 283 integer d_tr1(9),d_tr2(9),d_rr1(3),d_rr2(3),d_ir1(3),d_ir2(3) 284 integer k_f1_offset,k_v2_offset,k_d1_offset(3) 285 integer k_t1_offset,k_t2_offset 286 integer k_tr1_offset(3),k_tr2_offset(3) 287 integer size_tr1(3),size_tr2(3) 288 integer sym_abelian_axis 289 external sym_abelian_axis 290 double precision rr1,rr2,residual 291 double precision omega 292 double precision cpu, wall 293 double precision ddotfile 294 external ddotfile 295 logical nodezero 296 character*4 irrepname 297 character*255 filename 298 character*4 ir1filename(3) ! File name stub 299 data ir1filename/'ir1x','ir1y','ir1z'/ 300 character*4 ir2filename(3) ! File name stub 301 data ir2filename/'ir2x','ir2y','ir1z'/ 302 character*5 rr1filename(3) ! File name stub 303 data rr1filename/'rr1x ','rr1y ','rr1z '/ 304 character*5 rr2filename(3) ! File name stub 305 data rr2filename/'rr2x ','rr2y ','rr2z '/ 306c 307 nodezero=(ga_nodeid().eq.0) 308c 309 call tce_diis_init() 310 do iter=1,maxiter 311 cpu=-util_cpusec() 312 wall=-util_wallsec() 313 if (nodezero.and.(iter.eq.1)) 314 & write(LuOut,9400) "CCSD-IR (real component) - squared" 315c 316 call tce_filename(ir1filename(axis),filename) 317 call createfile(filename,d_ir1(axis),size_tr1(axis)) 318 call tce_zero(d_ir1(axis),size_tr1(axis)) 319c 320 call tce_filename(ir2filename(axis),filename) 321 call createfile(filename,d_ir2(axis),size_tr2(axis)) 322 call tce_zero(d_ir2(axis),size_tr2(axis)) 323c 324 call ccsd_o1(d_ir1(axis),d_d1(axis),d_t1,d_t2, 325 1 k_tr1_offset(axis),k_d1_offset(axis), 326 2 k_t1_offset,k_t2_offset) 327 call eomccsd_x1_old(d_f1,d_ir1(axis),d_t1,d_t2,d_v2, 328 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis),k_f1_offset, 329 2 k_tr1_offset(axis), 330 3 k_t1_offset,k_t2_offset,k_v2_offset, 331 4 k_tr1_offset(axis),k_tr2_offset(axis)) 332 call reconcilefile(d_ir1(axis),size_tr1(axis)) 333c 334 call ccsd_o2(d_ir2(axis),d_d1(axis),d_t1,d_t2, 335 1 k_tr2_offset(axis),k_d1_offset(axis), 336 2 k_t1_offset,k_t2_offset) 337 call eomccsd_x2_old(d_f1,d_ir2(axis),d_t1,d_t2,d_v2, 338 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis), 339 2 k_f1_offset,k_tr2_offset(axis), 340 3 k_t1_offset,k_t2_offset,k_v2_offset, 341 4 k_tr1_offset(axis),k_tr2_offset(axis), 342 5 size_tr1(axis),size_tr2(axis)) 343 call reconcilefile(d_ir2(axis),size_tr2(axis)) 344c 345 call tce_filename(rr1filename(axis),filename) 346 call createfile(filename,d_rr1(axis),size_tr1(axis)) 347 call tce_zero(d_rr1(axis),size_tr1(axis)) 348c 349 call tce_filename(rr2filename(axis),filename) 350 call createfile(filename,d_rr2(axis),size_tr2(axis)) 351 call tce_zero(d_rr2(axis),size_tr2(axis)) 352c 353 call eomccsd_x1_old(d_f1,d_rr1(axis),d_t1,d_t2,d_v2, 354 1 d_ir1(axis),d_ir2(axis),k_f1_offset, 355 2 k_tr1_offset(axis), 356 3 k_t1_offset,k_t2_offset,k_v2_offset, 357 4 k_tr1_offset(axis),k_tr2_offset(axis)) 358 call reconcilefile(d_rr1(axis),size_tr1(axis)) 359c 360 call eomccsd_x2_old(d_f1,d_rr2(axis),d_t1,d_t2,d_v2, 361 1 d_ir1(axis),d_ir2(axis), 362 2 k_f1_offset,k_tr2_offset(axis), 363 3 k_t1_offset,k_t2_offset,k_v2_offset, 364 4 k_tr1_offset(axis),k_tr2_offset(axis), 365 5 size_tr1(axis),size_tr2(axis)) 366 call reconcilefile(d_rr2(axis),size_tr2(axis)) 367c 368 call deletefile(d_ir1(axis)) 369 call deletefile(d_ir2(axis)) 370c 371 call daxpyfile(1,omega*omega,d_tr1(axis+dynaxis), 372 1 d_rr1(axis),size_tr1(axis)) 373 call reconcilefile(d_rr1(axis),size_tr1(axis)) 374 call tce_residual_tr1(d_rr1(axis),k_tr1_offset(axis),rr1) 375c 376 call daxpyfile(1,omega*omega,d_tr2(axis+dynaxis), 377 1 d_rr2(axis),size_tr2(axis)) 378 call reconcilefile(d_rr2(axis),size_tr2(axis)) 379 call tce_residual_tr2(d_rr2(axis),k_tr2_offset(axis),rr2) 380c 381 residual = max(rr1,rr2) 382 cpu=cpu+util_cpusec() 383 wall=wall+util_wallsec() 384 if (nodezero) write(LuOut,9420) iter,residual,cpu,wall 385 if (residual .lt. thresh) then 386 if (nodezero) write(LuOut,9410) 387 if (ampnorms) then 388 call tce_residual_tr1(d_tr1(axis+dynaxis), 389 1 k_tr1_offset(axis),rr1) 390 call tce_residual_tr2(d_tr2(axis+dynaxis), 391 1 k_tr2_offset(axis),rr2) 392 if (nodezero) then 393 write(LuOut,9082) "T(1) singles",rr1 394 write(LuOut,9082) "T(1) doubles",rr2 395 endif 396 endif 397 call deletefile(d_rr2(axis)) 398 call deletefile(d_rr1(axis)) 399 call tce_diis_tidy() 400 if (save_tr(1)) then 401 if(nodezero) then 402 write(LuOut,*) 'Saving T1(1) now...' 403 endif 404 call tr1_restart_save(d_tr1(axis+dynaxis), 405 1 k_tr1_offset(axis),size_tr1(axis), 406 2 axis+dynaxis,handle_tr1(axis),irrep_x) 407 endif 408 if (save_tr(2)) then 409 if(nodezero) then 410 write(LuOut,*) 'Saving T2(1) now...' 411 endif 412 call tr2_restart_save(d_tr2(axis+dynaxis), 413 1 k_tr2_offset(axis),size_tr2(axis), 414 2 axis+dynaxis,handle_tr2(axis),irrep_x) 415 endif 416 return 417 endif 418 if (save_tr(1).and.(mod(iter,save_interval).eq.0)) then 419 if(nodezero) then 420 write(LuOut,*) 'Saving T1(1) now...' 421 endif 422 call tr1_restart_save(d_tr1(axis+dynaxis), 423 1 k_tr1_offset(axis),size_tr1(axis), 424 2 axis+dynaxis,handle_tr1(axis),irrep_x) 425 endif 426 if (save_tr(2).and.(mod(iter,save_interval).eq.0)) then 427 if(nodezero) then 428 write(LuOut,*) 'Saving T2(1) now...' 429 endif 430 call tr2_restart_save(d_tr2(axis+dynaxis), 431 1 k_tr2_offset(axis),size_tr2(axis), 432 2 axis+dynaxis,handle_tr2(axis),irrep_x) 433 endif 434 call tce_diis3(.false.,iter,.true.,.true.,.false.,.false., 435 1 d_rr1(axis),d_tr1(axis+dynaxis),k_tr1_offset(axis), 436 2 size_tr1(axis), 437 3 d_rr2(axis),d_tr2(axis+dynaxis),k_tr2_offset(axis), 438 4 size_tr2(axis), 439 5 dummy,dummy,dummy,dummy, 440 6 dummy,dummy,dummy,dummy,omega,2) 441 call deletefile(d_rr2(axis)) 442 call deletefile(d_rr1(axis)) 443 if (nodezero) call util_flush(LuOut) 444 enddo ! iter loop 445 call errquit('ccsd_ir_real_sq_it: maxiter exceeded',iter,CALC_ERR) 446c 447 9020 format(1x,'Cpu & wall time / sec',2f15.1) 448 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 449 9100 format(1x,i4,2f18.13,2f8.1) 450 9120 format(1x,A) 451 9121 format(/,1x,A) 452 9122 format(1x,A,i4) 453 9400 format(/,1x,A,' iterations',/, 454 1 1x,'---------------------------------------------',/ 455 2 1x,'Iter Residuum Cpu Wall',/ 456 3 1x,'---------------------------------------------') 457 9410 format( 458 1 1x,'---------------------------------------------',/ 459 2 1x,'Iterations converged') 460 9420 format(1x,i4,f25.13,2f8.1) 461 9431 format(/,1x,'Frequency = ',f15.7,' / au') 462 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 463 return 464 end 465 466 467 subroutine ccsd_ir_imag_sq_it(axis,dynaxis,omega, 468 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 469 2 k_f1_offset,k_v2_offset,k_d1_offset, 470 3 k_t1_offset,k_t2_offset,k_tr1_offset,k_tr2_offset, 471 4 size_tr1,size_tr2) 472 implicit none 473#include "global.fh" 474#include "mafdecls.fh" 475#include "util.fh" 476#include "errquit.fh" 477#include "stdio.fh" 478#include "tce.fh" 479#include "tce_main.fh" 480#include "tce_diis.fh" 481#include "tce_prop.fh" 482#include "tce_restart.fh" 483c 484 integer i,j,dummy,axis 485 integer omegacount,omegasign 486 integer dynfreq,dynaxis 487 integer axisA,axisB 488 integer irrep_g 489 parameter (irrep_g=0) 490 integer d_f1,d_v2,d_d1(3),d_t1,d_t2 491 integer d_tr1(9),d_tr2(9),d_rr1(3),d_rr2(3),d_ir1(3),d_ir2(3) 492 integer k_f1_offset,k_v2_offset,k_d1_offset(3) 493 integer k_t1_offset,k_t2_offset 494 integer k_tr1_offset(3),k_tr2_offset(3) 495 integer size_tr1(3),size_tr2(3) 496 integer sym_abelian_axis 497 external sym_abelian_axis 498 double precision rr1,rr2,residual 499 double precision omega 500 double precision cpu, wall 501 double precision ddotfile 502 external ddotfile 503 logical nodezero 504 character*4 irrepname 505 character*255 filename 506 character*4 ir1filename(3) ! File name stub 507 data ir1filename/'ir1x','ir1y','ir1z'/ 508 character*4 ir2filename(3) ! File name stub 509 data ir2filename/'ir2x','ir2y','ir1z'/ 510 character*5 rr1filename(3) ! File name stub 511 data rr1filename/'rr1x ','rr1y ','rr1z '/ 512 character*5 rr2filename(3) ! File name stub 513 data rr2filename/'rr2x ','rr2y ','rr2z '/ 514c 515 nodezero=(ga_nodeid().eq.0) 516c 517 call tce_diis_init() 518 do iter=1,maxiter 519 cpu=-util_cpusec() 520 wall=-util_wallsec() 521 if (nodezero.and.(iter.eq.1)) 522 & write(LuOut,9400) "CCSD-IR (imaginary component) - squared" 523c 524 call tce_filename(ir1filename(axis),filename) 525 call createfile(filename,d_ir1(axis),size_tr1(axis)) 526 call tce_zero(d_ir1(axis),size_tr1(axis)) 527c 528 call tce_filename(ir2filename(axis),filename) 529 call createfile(filename,d_ir2(axis),size_tr2(axis)) 530 call tce_zero(d_ir2(axis),size_tr2(axis)) 531c 532 call eomccsd_x1_old(d_f1,d_ir1(axis),d_t1,d_t2,d_v2, 533 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis),k_f1_offset, 534 2 k_tr1_offset(axis), 535 3 k_t1_offset,k_t2_offset,k_v2_offset, 536 4 k_tr1_offset(axis),k_tr2_offset(axis)) 537 call reconcilefile(d_ir1(axis),size_tr1(axis)) 538c 539 call eomccsd_x2_old(d_f1,d_ir2(axis),d_t1,d_t2,d_v2, 540 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis), 541 2 k_f1_offset,k_tr2_offset(axis), 542 3 k_t1_offset,k_t2_offset,k_v2_offset, 543 4 k_tr1_offset(axis),k_tr2_offset(axis), 544 5 size_tr1(axis),size_tr2(axis)) 545 call reconcilefile(d_ir2(axis),size_tr2(axis)) 546c 547 call tce_filename(rr1filename(axis),filename) 548 call createfile(filename,d_rr1(axis),size_tr1(axis)) 549 call tce_zero(d_rr1(axis),size_tr1(axis)) 550c 551 call tce_filename(rr2filename(axis),filename) 552 call createfile(filename,d_rr2(axis),size_tr2(axis)) 553 call tce_zero(d_rr2(axis),size_tr2(axis)) 554c 555 call eomccsd_x1_old(d_f1,d_rr1(axis),d_t1,d_t2,d_v2, 556 1 d_ir1(axis),d_ir2(axis),k_f1_offset, 557 2 k_tr1_offset(axis),k_t1_offset,k_t2_offset,k_v2_offset, 558 4 k_tr1_offset(axis),k_tr2_offset(axis)) 559 call reconcilefile(d_rr1(axis),size_tr1(axis)) 560c 561 call eomccsd_x2_old(d_f1,d_rr2(axis),d_t1,d_t2,d_v2, 562 1 d_ir1(axis),d_ir2(axis),k_f1_offset, 563 2 k_tr2_offset(axis),k_t1_offset,k_t2_offset,k_v2_offset, 564 4 k_tr1_offset(axis),k_tr2_offset(axis), 565 5 size_tr1(axis),size_tr2(axis)) 566 call reconcilefile(d_rr2(axis),size_tr2(axis)) 567c 568 call tce_zero(d_ir1(axis),size_tr1(axis)) 569 call ccsd_o1(d_ir1(axis),d_d1(axis),d_t1,d_t2, 570 1 k_tr1_offset(axis),k_d1_offset(axis), 571 2 k_t1_offset,k_t2_offset) 572 call daxpyfile(1,omega,d_ir1(axis),d_rr1(axis),size_tr1(axis)) 573 call reconcilefile(d_rr1(axis),size_tr1(axis)) 574c 575 call tce_zero(d_ir2(axis),size_tr2(axis)) 576 call ccsd_o2(d_ir2(axis),d_d1(axis),d_t1,d_t2, 577 1 k_tr2_offset(axis),k_d1_offset(axis), 578 2 k_t1_offset,k_t2_offset) 579 call daxpyfile(1,omega,d_ir2(axis),d_rr2(axis),size_tr2(axis)) 580 call reconcilefile(d_rr2(axis),size_tr2(axis)) 581c 582 call deletefile(d_ir1(axis)) 583 call deletefile(d_ir2(axis)) 584c 585 call daxpyfile(1,omega*omega,d_tr1(axis+dynaxis), 586 1 d_rr1(axis),size_tr1(axis)) 587 call reconcilefile(d_rr1(axis),size_tr1(axis)) 588 call tce_residual_tr1(d_rr1(axis),k_tr1_offset(axis),rr1) 589c 590 call daxpyfile(1,omega*omega,d_tr2(axis+dynaxis), 591 1 d_rr2(axis),size_tr2(axis)) 592 call reconcilefile(d_rr2(axis),size_tr2(axis)) 593 call tce_residual_tr2(d_rr2(axis),k_tr2_offset(axis),rr2) 594c 595 residual = max(rr1,rr2) 596 cpu=cpu+util_cpusec() 597 wall=wall+util_wallsec() 598 if (nodezero) write(LuOut,9420) iter,residual,cpu,wall 599 if (residual .lt. thresh) then 600 if (nodezero) write(LuOut,9410) 601 if (ampnorms) then 602 call tce_residual_tr1(d_tr1(axis+dynaxis), 603 1 k_tr1_offset(axis),rr1) 604 call tce_residual_tr2(d_tr2(axis+dynaxis), 605 1 k_tr2_offset(axis),rr2) 606 if (nodezero) then 607 write(LuOut,9082) "T(1) singles",rr1 608 write(LuOut,9082) "T(1) doubles",rr2 609 endif 610 endif 611 call deletefile(d_rr2(axis)) 612 call deletefile(d_rr1(axis)) 613 call tce_diis_tidy() 614 if (save_tr(1)) then 615 if(nodezero) then 616 write(LuOut,*) 'Saving T1(1) now...' 617 endif 618 call tr1_restart_save(d_tr1(axis+dynaxis), 619 1 k_tr1_offset(axis),size_tr1(axis), 620 2 axis+dynaxis,handle_tr1(axis),irrep_x) 621 endif 622 if (save_tr(2)) then 623 if(nodezero) then 624 write(LuOut,*) 'Saving T2(1) now...' 625 endif 626 call tr2_restart_save(d_tr2(axis+dynaxis), 627 1 k_tr2_offset(axis),size_tr2(axis), 628 2 axis+dynaxis,handle_tr2(axis),irrep_x) 629 endif 630 return 631 endif 632 if (save_tr(1).and.(mod(iter,save_interval).eq.0)) then 633 if(nodezero) then 634 write(LuOut,*) 'Saving T1(1) now...' 635 endif 636 call tr1_restart_save(d_tr1(axis+dynaxis), 637 1 k_tr1_offset(axis),size_tr1(axis), 638 2 axis+dynaxis,handle_tr1(axis),irrep_x) 639 endif 640 if (save_tr(2).and.(mod(iter,save_interval).eq.0)) then 641 if(nodezero) then 642 write(LuOut,*) 'Saving T2(1) now...' 643 endif 644 call tr2_restart_save(d_tr2(axis+dynaxis), 645 1 k_tr2_offset(axis),size_tr2(axis), 646 2 axis+dynaxis,handle_tr2(axis),irrep_x) 647 endif 648 call tce_diis3(.false.,iter,.true.,.true.,.false.,.false., 649 1 d_rr1(axis),d_tr1(axis+dynaxis),k_tr1_offset(axis), 650 2 size_tr1(axis), 651 3 d_rr2(axis),d_tr2(axis+dynaxis),k_tr2_offset(axis), 652 4 size_tr2(axis), 653 5 dummy,dummy,dummy,dummy, 654 6 dummy,dummy,dummy,dummy,omega,2) 655 call deletefile(d_rr2(axis)) 656 call deletefile(d_rr1(axis)) 657 if (nodezero) call util_flush(LuOut) 658 enddo ! iter loop 659 call errquit('ccsd_ir_imag_sq_it: maxiter exceeded',iter,CALC_ERR) 660c 661 9020 format(1x,'Cpu & wall time / sec',2f15.1) 662 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 663 9100 format(1x,i4,2f18.13,2f8.1) 664 9120 format(1x,A) 665 9121 format(/,1x,A) 666 9122 format(1x,A,i4) 667 9400 format(/,1x,A,' iterations',/, 668 1 1x,'---------------------------------------------',/ 669 2 1x,'Iter Residuum Cpu Wall',/ 670 3 1x,'---------------------------------------------') 671 9410 format( 672 1 1x,'---------------------------------------------',/ 673 2 1x,'Iterations converged') 674 9420 format(1x,i4,f25.13,2f8.1) 675 9431 format(/,1x,'Frequency = ',f15.7,' / au') 676 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 677 return 678 end 679 680 681 682 subroutine ccsd_ir_real_it(axis,dynaxis,omega, 683 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 684 2 k_f1_offset,k_v2_offset,k_d1_offset, 685 3 k_t1_offset,k_t2_offset, 686 4 k_tr1_offset,k_tr2_offset, 687 5 size_tr1,size_tr2) 688 implicit none 689#include "global.fh" 690#include "mafdecls.fh" 691#include "util.fh" 692#include "errquit.fh" 693#include "stdio.fh" 694#include "tce.fh" 695#include "tce_main.fh" 696#include "tce_diis.fh" 697#include "tce_prop.fh" 698#include "tce_restart.fh" 699c 700 integer i,j,dummy,axis 701 integer omegacount 702 integer omegasign 703 integer dynfreq 704 integer dynaxis 705 integer irrep_g 706 parameter (irrep_g=0) 707 integer d_f1,d_v2,d_d1(3),d_t1,d_t2 708 integer d_tr1(9),d_tr2(9),d_rr1(3),d_rr2(3),d_ir1(3),d_ir2(3) 709 integer k_f1_offset,k_v2_offset,k_d1_offset(3) 710 integer k_t1_offset,k_t2_offset 711 integer k_tr1_offset(3),k_tr2_offset(3) 712 integer size_tr1(3),size_tr2(3) 713 integer sym_abelian_axis 714 external sym_abelian_axis 715 double precision rr1,rr2,residual 716 double precision omega 717 double precision cpu, wall 718 double precision ddotfile 719 external ddotfile 720 logical nodezero 721 character*4 irrepname 722 character*255 filename 723 character*4 ir1filename(3) ! File name stub 724 data ir1filename/'ir1x','ir1y','ir1z'/ 725 character*4 ir2filename(3) ! File name stub 726 data ir2filename/'ir2x','ir2y','ir1z'/ 727 character*5 rr1filename(3) ! File name stub 728 data rr1filename/'rr1x ','rr1y ','rr1z '/ 729 character*5 rr2filename(3) ! File name stub 730 data rr2filename/'rr2x ','rr2y ','rr2z '/ 731c 732 nodezero=(ga_nodeid().eq.0) 733c 734 call tce_diis_init() 735 do iter=1,maxiter 736 cpu=-util_cpusec() 737 wall=-util_wallsec() 738 if (nodezero.and.(iter.eq.1)) 739 & write(LuOut,9400) "CCSD-IR (real component) - linear" 740c 741 call tce_filename(rr1filename(axis),filename) 742 call createfile(filename,d_rr1(axis),size_tr1(axis)) 743 call tce_zero(d_rr1(axis),size_tr1(axis)) 744c 745 call tce_filename(rr2filename(axis),filename) 746 call createfile(filename,d_rr2(axis),size_tr2(axis)) 747 call tce_zero(d_rr2(axis),size_tr2(axis)) 748c 749 call ccsd_o1(d_rr1(axis),d_d1(axis),d_t1,d_t2, 750 1 k_tr1_offset(axis),k_d1_offset(axis), 751 2 k_t1_offset,k_t2_offset) 752 call daxfile(1,(-1.0d0)*omega,d_tr1(axis), 753 1 d_rr1(axis),size_tr1(axis)) 754 call eomccsd_x1_old(d_f1,d_rr1(axis),d_t1,d_t2,d_v2, 755 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis), 756 2 k_f1_offset,k_tr1_offset(axis), 757 3 k_t1_offset,k_t2_offset,k_v2_offset, 758 4 k_tr1_offset(axis),k_tr2_offset(axis)) 759 call reconcilefile(d_rr1(axis),size_tr1(axis)) 760 call tce_residual_tr1(d_rr1(axis),k_tr1_offset(axis),rr1) 761c 762 call ccsd_o2(d_rr2(axis),d_d1(axis),d_t1,d_t2, 763 1 k_tr2_offset(axis),k_d1_offset(axis), 764 2 k_t1_offset,k_t2_offset) 765 call daxfile(1,(-1.0d0)*omega,d_tr2(axis), 766 1 d_rr2(axis),size_tr2(axis)) 767 call eomccsd_x2_old(d_f1,d_rr2(axis),d_t1,d_t2,d_v2, 768 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis), 769 2 k_f1_offset,k_tr2_offset(axis), 770 3 k_t1_offset,k_t2_offset,k_v2_offset, 771 4 k_tr1_offset(axis),k_tr2_offset(axis), 772 5 size_tr1(axis),size_tr2(axis)) 773 call reconcilefile(d_rr2(axis),size_tr2(axis)) 774 call tce_residual_tr2(d_rr2(axis),k_tr2_offset(axis),rr2) 775c 776 residual = max(rr1,rr2) 777 cpu=cpu+util_cpusec() 778 wall=wall+util_wallsec() 779 if (nodezero) write(LuOut,9420) iter,residual,cpu,wall 780 if (residual .lt. thresh) then 781 if (nodezero) write(LuOut,9410) 782 if (ampnorms) then 783 call tce_residual_tr1(d_tr1(axis+dynaxis), 784 1 k_tr1_offset(axis),rr1) 785 call tce_residual_tr2(d_tr2(axis+dynaxis), 786 1 k_tr2_offset(axis),rr2) 787 if (nodezero) then 788 write(LuOut,9082) "T(1) singles",rr1 789 write(LuOut,9082) "T(1) doubles",rr2 790 endif 791 endif 792 call deletefile(d_rr2(axis)) 793 call deletefile(d_rr1(axis)) 794 call tce_diis_tidy() 795 if (save_tr(1)) then 796 if(nodezero) then 797 write(LuOut,*) 'Saving T1(1) now...' 798 endif 799 call tr1_restart_save(d_tr1(axis+dynaxis), 800 1 k_tr1_offset(axis),size_tr1(axis), 801 2 axis+dynaxis,handle_tr1(axis),irrep_x) 802 endif 803 if (save_tr(2)) then 804 if(nodezero) then 805 write(LuOut,*) 'Saving T2(1) now...' 806 endif 807 call tr2_restart_save(d_tr2(axis+dynaxis), 808 1 k_tr2_offset(axis),size_tr2(axis), 809 2 axis+dynaxis,handle_tr2(axis),irrep_x) 810 endif 811 return 812 endif 813 if (save_tr(1).and.(mod(iter,save_interval).eq.0)) then 814 if(nodezero) then 815 write(LuOut,*) 'Saving T1(1) now...' 816 endif 817 call tr1_restart_save(d_tr1(axis+dynaxis), 818 1 k_tr1_offset(axis),size_tr1(axis), 819 2 axis+dynaxis,handle_tr1(axis),irrep_x) 820 endif 821 if (save_tr(2).and.(mod(iter,save_interval).eq.0)) then 822 if(nodezero) then 823 write(LuOut,*) 'Saving T2(1) now...' 824 endif 825 call tr2_restart_save(d_tr2(axis+dynaxis), 826 1 k_tr2_offset(axis),size_tr2(axis), 827 2 axis+dynaxis,handle_tr2(axis),irrep_x) 828 endif 829 call tce_diis3(.false.,iter,.true.,.true.,.false.,.false., 830 1 d_rr1(axis),d_tr1(axis+dynaxis),k_tr1_offset(axis), 831 2 size_tr1(axis), 832 3 d_rr2(axis),d_tr2(axis+dynaxis),k_tr2_offset(axis), 833 4 size_tr2(axis), 834 5 dummy,dummy,dummy,dummy, 835 6 dummy,dummy,dummy,dummy,0.0d0,1) 836 call deletefile(d_rr2(axis)) 837 call deletefile(d_rr1(axis)) 838 if (nodezero) call util_flush(LuOut) 839 enddo ! iter loop 840 call errquit('ccsd_ir_real_it: maxiter exceeded',iter,CALC_ERR) 841c 842 9020 format(1x,'Cpu & wall time / sec',2f15.1) 843 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 844 9100 format(1x,i4,2f18.13,2f8.1) 845 9120 format(1x,A) 846 9121 format(/,1x,A) 847 9122 format(1x,A,i4) 848 9400 format(/,1x,A,' iterations',/, 849 1 1x,'---------------------------------------------',/ 850 2 1x,'Iter Residuum Cpu Wall',/ 851 3 1x,'---------------------------------------------') 852 9410 format( 853 1 1x,'---------------------------------------------',/ 854 2 1x,'Iterations converged') 855 9420 format(1x,i4,f25.13,2f8.1) 856 9431 format(/,1x,'Frequency = ',f15.7,' / au') 857 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 858 return 859 end 860 861 862 subroutine ccsd_ir_imag_it(axis,dynaxis,omega, 863 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 864 2 k_f1_offset,k_v2_offset,k_d1_offset, 865 3 k_t1_offset,k_t2_offset, 866 4 k_tr1_offset,k_tr2_offset, 867 5 size_tr1,size_tr2) 868 implicit none 869#include "global.fh" 870#include "mafdecls.fh" 871#include "util.fh" 872#include "errquit.fh" 873#include "stdio.fh" 874#include "tce.fh" 875#include "tce_main.fh" 876#include "tce_diis.fh" 877#include "tce_prop.fh" 878#include "tce_restart.fh" 879c 880 integer i,j,dummy,axis 881 integer omegacount 882 integer omegasign 883 integer dynfreq 884 integer dynaxis 885 integer irrep_g 886 parameter (irrep_g=0) 887 integer d_f1,d_v2,d_d1(3),d_t1,d_t2,d_tr1(9),d_tr2(9) 888 integer d_rr1(3),d_rr2(3),d_ir1(3),d_ir2(3) 889 integer k_f1_offset,k_v2_offset,k_d1_offset(3) 890 integer k_t1_offset,k_t2_offset 891 integer k_tr1_offset(3),k_tr2_offset(3) 892 integer size_tr1(3),size_tr2(3) 893 integer sym_abelian_axis 894 external sym_abelian_axis 895 double precision rr1,rr2,residual 896 double precision omega 897 double precision cpu, wall 898 double precision ddotfile 899 external ddotfile 900 logical nodezero 901 character*4 irrepname 902 character*255 filename 903 character*4 ir1filename(3) ! File name stub 904 data ir1filename/'ir1x','ir1y','ir1z'/ 905 character*4 ir2filename(3) ! File name stub 906 data ir2filename/'ir2x','ir2y','ir1z'/ 907 character*5 rr1filename(3) ! File name stub 908 data rr1filename/'rr1x ','rr1y ','rr1z '/ 909 character*5 rr2filename(3) ! File name stub 910 data rr2filename/'rr2x ','rr2y ','rr2z '/ 911c 912 nodezero=(ga_nodeid().eq.0) 913c 914 call tce_diis_init() 915 do iter=1,maxiter 916 cpu=-util_cpusec() 917 wall=-util_wallsec() 918 if (nodezero.and.(iter.eq.1)) 919 & write(LuOut,9400) "CCSD-IR (imaginary component) - linear" 920c 921 call tce_filename(rr1filename(axis),filename) 922 call createfile(filename,d_rr1(axis),size_tr1(axis)) 923 call tce_zero(d_rr1(axis),size_tr1(axis)) 924c 925 call tce_filename(rr2filename(axis),filename) 926 call createfile(filename,d_rr2(axis),size_tr2(axis)) 927 call tce_zero(d_rr2(axis),size_tr2(axis)) 928c 929 call daxfile(1,(-1.0d0)*omega,d_tr1(axis+6), 930 1 d_rr1(axis),size_tr1(axis)) 931 call eomccsd_x1_old(d_f1,d_rr1(axis),d_t1,d_t2,d_v2, 932 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis), 933 2 k_f1_offset,k_tr1_offset(axis), 934 3 k_t1_offset,k_t2_offset,k_v2_offset, 935 4 k_tr1_offset(axis),k_tr2_offset(axis)) 936 call reconcilefile(d_rr1(axis),size_tr1(axis)) 937 call tce_residual_tr1(d_rr1(axis),k_tr1_offset(axis),rr1) 938c 939 call daxfile(1,(-1.0d0)*omega,d_tr2(axis+6), 940 1 d_rr2(axis),size_tr2(axis)) 941 call eomccsd_x2_old(d_f1,d_rr2(axis),d_t1,d_t2,d_v2, 942 1 d_tr1(axis+dynaxis),d_tr2(axis+dynaxis), 943 2 k_f1_offset,k_tr2_offset(axis), 944 3 k_t1_offset,k_t2_offset,k_v2_offset, 945 4 k_tr1_offset(axis),k_tr2_offset(axis), 946 5 size_tr1(axis),size_tr2(axis)) 947 call reconcilefile(d_rr2(axis),size_tr2(axis)) 948 call tce_residual_tr2(d_rr2(axis),k_tr2_offset(axis),rr2) 949c 950 residual = max(rr1,rr2) 951 cpu=cpu+util_cpusec() 952 wall=wall+util_wallsec() 953 if (nodezero) write(LuOut,9420) iter,residual,cpu,wall 954 if (residual .lt. thresh) then 955 if (nodezero) write(LuOut,9410) 956 if (ampnorms) then 957 call tce_residual_tr1(d_tr1(axis+dynaxis), 958 1 k_tr1_offset(axis),rr1) 959 call tce_residual_tr2(d_tr2(axis+dynaxis), 960 1 k_tr2_offset(axis),rr2) 961 if (nodezero) then 962 write(LuOut,9082) "T(1) singles",rr1 963 write(LuOut,9082) "T(1) doubles",rr2 964 endif 965 endif 966 call deletefile(d_rr2(axis)) 967 call deletefile(d_rr1(axis)) 968 call tce_diis_tidy() 969 if (save_tr(1)) then 970 if(nodezero) then 971 write(LuOut,*) 'Saving T1(1) now...' 972 endif 973 call tr1_restart_save(d_tr1(axis+dynaxis), 974 1 k_tr1_offset(axis),size_tr1(axis), 975 2 axis+dynaxis,handle_tr1(axis),irrep_x) 976 endif 977 if (save_tr(2)) then 978 if(nodezero) then 979 write(LuOut,*) 'Saving T2(1) now...' 980 endif 981 call tr2_restart_save(d_tr2(axis+dynaxis), 982 1 k_tr2_offset(axis),size_tr2(axis), 983 2 axis+dynaxis,handle_tr2(axis),irrep_x) 984 endif 985 return 986 endif 987 if (save_tr(1).and.(mod(iter,save_interval).eq.0)) then 988 if(nodezero) then 989 write(LuOut,*) 'Saving T1(1) now...' 990 endif 991 call tr1_restart_save(d_tr1(axis+dynaxis), 992 1 k_tr1_offset(axis),size_tr1(axis), 993 2 axis+dynaxis,handle_tr1(axis),irrep_x) 994 endif 995 if (save_tr(2).and.(mod(iter,save_interval).eq.0)) then 996 if(nodezero) then 997 write(LuOut,*) 'Saving T2(1) now...' 998 endif 999 call tr2_restart_save(d_tr2(axis+dynaxis), 1000 1 k_tr2_offset(axis),size_tr2(axis), 1001 2 axis+dynaxis,handle_tr2(axis),irrep_x) 1002 endif 1003 call tce_diis3(.false.,iter,.true.,.true.,.false.,.false., 1004 1 d_rr1(axis),d_tr1(axis+dynaxis),k_tr1_offset(axis), 1005 2 size_tr1(axis), 1006 3 d_rr2(axis),d_tr2(axis+dynaxis),k_tr2_offset(axis), 1007 4 size_tr2(axis), 1008 5 dummy,dummy,dummy,dummy, 1009 6 dummy,dummy,dummy,dummy,0.0d0,1) 1010 call deletefile(d_rr2(axis)) 1011 call deletefile(d_rr1(axis)) 1012 if (nodezero) call util_flush(LuOut) 1013 enddo ! iter loop 1014 call errquit('ccsd_ir_imag_it: maxiter exceeded',iter,CALC_ERR) 1015c 1016 9020 format(1x,'Cpu & wall time / sec',2f15.1) 1017 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 1018 9100 format(1x,i4,2f18.13,2f8.1) 1019 9120 format(1x,A) 1020 9121 format(/,1x,A) 1021 9122 format(1x,A,i4) 1022 9400 format(/,1x,A,' iterations',/, 1023 1 1x,'---------------------------------------------',/ 1024 2 1x,'Iter Residuum Cpu Wall',/ 1025 3 1x,'---------------------------------------------') 1026 9410 format( 1027 1 1x,'---------------------------------------------',/ 1028 2 1x,'Iterations converged') 1029 9420 format(1x,i4,f25.13,2f8.1) 1030 9431 format(/,1x,'Frequency = ',f15.7,' / au') 1031 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 1032 return 1033 end 1034 1035 1036 1037 subroutine ccsd_ir_coupled_it(guess,axis,omega, 1038 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 1039 2 k_f1_offset,k_v2_offset,k_d1_offset, 1040 3 k_t1_offset,k_t2_offset, 1041 4 k_tr1_offset,k_tr2_offset, 1042 5 size_tr1,size_tr2) 1043 implicit none 1044#include "global.fh" 1045#include "mafdecls.fh" 1046#include "util.fh" 1047#include "errquit.fh" 1048#include "stdio.fh" 1049#include "tce.fh" 1050#include "tce_main.fh" 1051#include "tce_diis.fh" 1052#include "tce_prop.fh" 1053#include "tce_restart.fh" 1054c 1055 integer i,j,dummy,axis,dynaxis 1056 integer omegacount 1057 integer omegasign 1058 integer irrep_g 1059 parameter (irrep_g=0) 1060 integer d_f1,d_v2,d_d1(3),d_t1,d_t2 1061 integer d_tr1(9),d_tr2(9) 1062 integer d_rr1(3),d_rr2(3),d_ir1(3),d_ir2(3) 1063 integer k_f1_offset,k_v2_offset,k_d1_offset(3) 1064 integer k_t1_offset,k_t2_offset 1065 integer k_tr1_offset(3),k_tr2_offset(3) 1066 integer size_tr1(3),size_tr2(3) 1067 integer sym_abelian_axis 1068 integer zero,one 1069 parameter(zero = 0) 1070 parameter(one = 1) 1071 external sym_abelian_axis 1072 double precision rr1,rr2,rr3,rr4,rres,ires,residual 1073 double precision omega 1074 double precision cpu, wall 1075 double precision ddotfile 1076 external ddotfile 1077 logical nodezero,guess 1078 character*4 irrepname 1079 character*255 filename 1080 character*5 rr1filename(3) ! File name stub 1081 data rr1filename/'rr1x ','rr1y ','rr1z '/ 1082 character*5 rr2filename(3) ! File name stub 1083 data rr2filename/'rr2x ','rr2y ','rr2z '/ 1084 character*4 ir1filename(3) ! File name stub 1085 data ir1filename/'ir1x','ir1y','ir1z'/ 1086 character*4 ir2filename(3) ! File name stub 1087 data ir2filename/'ir2x','ir2y','ir1z'/ 1088c 1089 nodezero=(ga_nodeid().eq.0) 1090c 1091 call tce_diis_init() 1092 do iter=1,maxiter 1093 cpu=-util_cpusec() 1094 wall=-util_wallsec() 1095c 1096c REAL COMPONENT 1097c 1098 if (guess.and.(iter.eq.1).and.guess_ir_real) then 1099 call ccsd_lr_guess(d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 1100 1 k_f1_offset,k_v2_offset,k_d1_offset, 1101 2 k_t1_offset,k_t2_offset,k_tr1_offset,k_tr2_offset, 1102 3 size_tr1,size_tr2,axis,6) 1103 endif 1104c 1105 if (guess.and.(iter.eq.1).and.guess_ir_imag) then 1106 call ccsd_ir_guess(d_d1,d_t1,d_t2,d_tr1,d_tr2, 1107 1 k_d1_offset,k_t1_offset,k_t2_offset, 1108 2 k_tr1_offset,k_tr2_offset, 1109 4 size_tr1,size_tr2,one,axis,omega) 1110 endif 1111c 1112 call tce_filename(rr1filename(axis),filename) 1113 call createfile(filename,d_rr1(axis),size_tr1(axis)) 1114 call tce_zero(d_rr1(axis),size_tr1(axis)) 1115c 1116 call tce_filename(rr2filename(axis),filename) 1117 call createfile(filename,d_rr2(axis),size_tr2(axis)) 1118 call tce_zero(d_rr2(axis),size_tr2(axis)) 1119c 1120 call daxfile(1,(-1.0d0)*omega,d_tr1(axis+0), 1121 1 d_rr1(axis),size_tr1(axis)) 1122 call ccsd_o1(d_rr1(axis),d_d1(axis),d_t1,d_t2, 1123 1 k_tr1_offset(axis),k_d1_offset(axis), 1124 2 k_t1_offset,k_t2_offset) 1125 call eomccsd_x1_old(d_f1,d_rr1(axis),d_t1,d_t2,d_v2, 1126 1 d_tr1(axis+6),d_tr2(axis+6), 1127 2 k_f1_offset,k_tr1_offset(axis), 1128 3 k_t1_offset,k_t2_offset,k_v2_offset, 1129 4 k_tr1_offset(axis),k_tr2_offset(axis)) 1130 call reconcilefile(d_rr1(axis),size_tr1(axis)) 1131 call tce_residual_tr1(d_rr1(axis),k_tr1_offset(axis),rr1) 1132c 1133 call daxfile(1,(-1.0d0)*omega,d_tr2(axis+0), 1134 1 d_rr2(axis),size_tr2(axis)) 1135 call ccsd_o2(d_rr2(axis),d_d1(axis),d_t1,d_t2, 1136 1 k_tr2_offset(axis),k_d1_offset(axis), 1137 2 k_t1_offset,k_t2_offset) 1138 call eomccsd_x2_old(d_f1,d_rr2(axis),d_t1,d_t2,d_v2, 1139 1 d_tr1(axis+6),d_tr2(axis+6), 1140 2 k_f1_offset,k_tr2_offset(axis), 1141 3 k_t1_offset,k_t2_offset,k_v2_offset, 1142 4 k_tr1_offset(axis),k_tr2_offset(axis), 1143 5 size_tr1(axis),size_tr2(axis)) 1144 call reconcilefile(d_rr2(axis),size_tr2(axis)) 1145 call tce_residual_tr2(d_rr2(axis),k_tr2_offset(axis),rr2) 1146c 1147c IMAGINARY COMPONENT 1148c 1149 call tce_filename(ir1filename(axis),filename) 1150 call createfile(filename,d_ir1(axis),size_tr1(axis)) 1151 call tce_zero(d_ir1(axis),size_tr1(axis)) 1152c 1153 call tce_filename(ir2filename(axis),filename) 1154 call createfile(filename,d_ir2(axis),size_tr2(axis)) 1155 call tce_zero(d_ir2(axis),size_tr2(axis)) 1156c 1157 call daxfile(1,(-1.0d0)*omega,d_tr1(axis+6), 1158 1 d_ir1(axis),size_tr1(axis)) 1159 call eomccsd_x1_old(d_f1,d_ir1(axis),d_t1,d_t2,d_v2, 1160 1 d_tr1(axis+0),d_tr2(axis+0), 1161 2 k_f1_offset,k_tr1_offset(axis), 1162 3 k_t1_offset,k_t2_offset,k_v2_offset, 1163 4 k_tr1_offset(axis),k_tr2_offset(axis)) 1164 call reconcilefile(d_ir1(axis),size_tr1(axis)) 1165 call tce_residual_tr1(d_ir1(axis),k_tr1_offset(axis),rr3) 1166c 1167 call daxfile(1,(-1.0d0)*omega,d_tr2(axis+6), 1168 1 d_ir2(axis),size_tr2(axis)) 1169 call eomccsd_x2_old(d_f1,d_ir2(axis),d_t1,d_t2,d_v2, 1170 1 d_tr1(axis+0),d_tr2(axis+0), 1171 2 k_f1_offset,k_tr2_offset(axis), 1172 3 k_t1_offset,k_t2_offset,k_v2_offset, 1173 4 k_tr1_offset(axis),k_tr2_offset(axis), 1174 5 size_tr1(axis),size_tr2(axis)) 1175 call reconcilefile(d_ir2(axis),size_tr2(axis)) 1176 call tce_residual_tr2(d_ir2(axis),k_tr2_offset(axis),rr4) 1177c 1178 if (nodezero.and.(iter.eq.1)) then 1179 write(LuOut,9050) "CCSD-IR (coupled iteration)" 1180 endif 1181c 1182 rres = max(rr1,rr2) 1183 ires = max(rr3,rr4) 1184 residual = max(rres,ires) 1185 cpu=cpu+util_cpusec() 1186 wall=wall+util_wallsec() 1187 if (nodezero) write(LuOut,9100) iter,rres,ires,cpu,wall 1188 if (residual .lt. thresh) then 1189 if (nodezero) write(LuOut,9060) 1190 call deletefile(d_rr2(axis)) 1191 call deletefile(d_rr1(axis)) 1192 call deletefile(d_ir2(axis)) 1193 call deletefile(d_ir1(axis)) 1194 call tce_diis_tidy() 1195 return 1196 endif 1197 call tce_diis3c(iter,.true.,.true.,.true.,.true., 1198 1 d_rr1(axis),d_tr1(axis+6),k_tr1_offset(axis),size_tr1(axis), 1199 2 d_rr2(axis),d_tr2(axis+6),k_tr2_offset(axis),size_tr2(axis), 1200 3 d_ir1(axis),d_tr1(axis+0),k_tr1_offset(axis),size_tr1(axis), 1201 4 d_ir2(axis),d_tr2(axis+0),k_tr2_offset(axis),size_tr2(axis), 1202 5 0.0d0) 1203 call deletefile(d_rr2(axis)) 1204 call deletefile(d_rr1(axis)) 1205 call deletefile(d_ir2(axis)) 1206 call deletefile(d_ir1(axis)) 1207 if (nodezero) call util_flush(LuOut) 1208 enddo ! iter loop 1209 call errquit('ccsd_ir_coupled_it: maxiter exceeded',iter,CALC_ERR) 1210c 1211 9020 format(1x,'Cpu & wall time / sec',2f15.1) 1212 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 1213 9120 format(1x,A) 1214 9121 format(/,1x,A) 1215 9122 format(1x,A,i4) 1216 9050 format(/,1x,A,' iterations',/, 1217 1 1x,'--------------------------------------------------------',/ 1218 2 1x,'Iter Residuum (real) Residuum (imag) Cpu Wall',/ 1219 3 1x,'--------------------------------------------------------') 1220 9060 format( 1221 1 1x,'--------------------------------------------------------',/ 1222 2 1x,'Iterations converged') 1223 9100 format(1x,i4,2f18.13,2f8.1) 1224 9431 format(/,1x,'Frequency = ',f15.7,' / au') 1225 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 1226 return 1227 end 1228 1229 1230 subroutine ccsd_ir_get_imag_from_real(axis,omega, 1231 1 d_f1,d_v2,d_d1,d_t1,d_t2,d_tr1,d_tr2, 1232 2 k_f1_offset,k_v2_offset,k_d1_offset, 1233 3 k_t1_offset,k_t2_offset,k_tr1_offset,k_tr2_offset, 1234 4 size_tr1,size_tr2) 1235 implicit none 1236#include "global.fh" 1237#include "mafdecls.fh" 1238#include "util.fh" 1239#include "errquit.fh" 1240#include "stdio.fh" 1241#include "tce.fh" 1242#include "tce_main.fh" 1243#include "tce_diis.fh" 1244#include "tce_prop.fh" 1245#include "tce_restart.fh" 1246c 1247 integer i,j,dummy,axis 1248 integer irrep_g 1249 parameter (irrep_g=0) 1250 integer d_f1,d_v2,d_d1(3),d_t1,d_t2 1251 integer d_tr1(9),d_tr2(9) 1252 integer k_f1_offset,k_v2_offset,k_d1_offset(3) 1253 integer k_t1_offset,k_t2_offset 1254 integer k_tr1_offset(3),k_tr2_offset(3) 1255 integer size_tr1(3),size_tr2(3) 1256 integer sym_abelian_axis 1257 external sym_abelian_axis 1258 double precision omega 1259 double precision cpu, wall 1260 double precision neg 1261 parameter (neg=-1.0d0) 1262 logical nodezero 1263 character*4 irrepname 1264 character*255 filename 1265c 1266 nodezero=(ga_nodeid().eq.0) 1267c 1268 cpu=-util_cpusec() 1269 wall=-util_wallsec() 1270c 1271 if (omega.eq.(0.0d0)) then 1272c 1273c x_I = 0 if w_I = 0 1274c 1275 call tce_zero(d_tr1(axis+0),size_tr1(axis)) 1276 call tce_zero(d_tr2(axis+0),size_tr2(axis)) 1277 call tce_zero(d_tr1(axis+3),size_tr1(axis)) 1278 call tce_zero(d_tr2(axis+3),size_tr2(axis)) 1279 return 1280 else 1281c 1282c x_I(+) = ( A * x_R - b ) / omega_I 1283c x_I(-) = -x_I(+) 1284c 1285 call tce_zero(d_tr1(axis+0),size_tr1(axis)) 1286 call tce_zero(d_tr2(axis+0),size_tr2(axis)) 1287c 1288 call ccsd_o1(d_tr1(axis),d_d1(axis),d_t1,d_t2, 1289 1 k_tr1_offset(axis),k_d1_offset(axis), 1290 2 k_t1_offset,k_t2_offset) 1291 call eomccsd_x1_old(d_f1,d_tr1(axis+0),d_t1,d_t2,d_v2, 1292 1 d_tr1(axis+6),d_tr2(axis+6),k_f1_offset, 1293 2 k_tr1_offset(axis), 1294 3 k_t1_offset,k_t2_offset,k_v2_offset, 1295 4 k_tr1_offset(axis),k_tr2_offset(axis)) 1296c 1297 call ccsd_o2(d_tr2(axis),d_d1(axis),d_t1,d_t2, 1298 1 k_tr2_offset(axis),k_d1_offset(axis), 1299 2 k_t1_offset,k_t2_offset) 1300 call eomccsd_x2_old(d_f1,d_tr2(axis+0),d_t1,d_t2,d_v2, 1301 1 d_tr1(axis+6),d_tr2(axis+6), 1302 2 k_f1_offset,k_tr2_offset(axis), 1303 3 k_t1_offset,k_t2_offset,k_v2_offset, 1304 4 k_tr1_offset(axis),k_tr2_offset(axis), 1305 5 size_tr1(axis),size_tr2(axis)) 1306c 1307 call reconcilefile(d_tr1(axis+0),size_tr1(axis)) 1308 call reconcilefile(d_tr2(axis+0),size_tr2(axis)) 1309c 1310 call dscalfile((1.0d0/omega),d_tr1(axis+0),size_tr1(axis)) 1311 call dscalfile((1.0d0/omega),d_tr2(axis+0),size_tr2(axis)) 1312c 1313c call daxfile(1,neg,d_tr1(axis+0),d_tr1(axis+3),size_tr1(axis)) 1314c call daxfile(1,neg,d_tr2(axis+0),d_tr2(axis+3),size_tr2(axis)) 1315c call copyfile(d_tr1(axis+0),d_tr1(axis+3),size_tr1(axis)) 1316c call copyfile(d_tr2(axis+0),d_tr2(axis+3),size_tr2(axis)) 1317 return 1318 endif 1319c 1320 9020 format(1x,'Cpu & wall time / sec',2f15.1) 1321 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 1322 9100 format(1x,i4,2f18.13,2f8.1) 1323 9120 format(1x,A) 1324 9121 format(/,1x,A) 1325 9122 format(1x,A,i4) 1326 9420 format(1x,i4,f25.13,2f8.1) 1327 9431 format(/,1x,'Frequency = ',f15.7,' / au') 1328 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 1329 end 1330 1331 1332 subroutine ccsd_ir_eval(omegacount,d_a0,d_f1,d_v2,d_d1, 1333 1 d_t1,d_t2,d_lambda1,d_lambda2,d_tr1,d_tr2, 1334 2 k_a0_offset,k_f1_offset,k_v2_offset,k_d1_offset, 1335 4 k_t1_offset,k_t2_offset,k_l1_offset,k_l2_offset, 1336 6 k_tr1_offset,k_tr2_offset) 1337 implicit none 1338#include "global.fh" 1339#include "mafdecls.fh" 1340#include "util.fh" 1341#include "errquit.fh" 1342#include "stdio.fh" 1343#include "tce.fh" 1344#include "tce_main.fh" 1345#include "tce_prop.fh" 1346c 1347 integer i,j,dummy,axis 1348 integer axisA 1349 integer axisB 1350 integer omegacount 1351 integer omegasign 1352 integer dynfreq 1353 integer dynaxis 1354 integer irrep_g 1355 parameter (irrep_g=0) 1356 integer d_a0,d_f1,d_v2,d_d1(3) 1357 integer d_t1,d_t2,d_lambda1,d_lambda2 1358 integer k_a0_offset,k_f1_offset,k_v2_offset,k_d1_offset(3) 1359 integer k_t1_offset,k_t2_offset,k_l1_offset,k_l2_offset 1360 integer d_tr1(9),d_tr2(9) 1361 integer k_tr1_offset(3),k_tr2_offset(3) 1362 integer sym_abelian_axis 1363 external sym_abelian_axis 1364 double precision omega 1365 double precision cpu, wall 1366 double precision alpha1,alpha2,alpha3,alpha4 1367 double precision au2ang ! Conversion factor from bohr to Angstrom 1368 double precision au2ang3 ! Conversion factor from bohr^3 to Angstrom^3 1369 double precision auXnm ! Conversion factor from a.u. (frequency) to nm (wavelength) 1370 double precision alpha(3,3) ! Dipole polarizability tensor 1371 double precision alphacopy(3,3) ! Dipole polarizability tensor copy 1372 double precision alphaiso ! Isotropic dipole polarizability 1373 double precision alphaani ! Anisotropic dipole polarizability 1374 double precision alphaevr(3) ! Dipole polarizability tensor eigenvalues (real) 1375 double precision alphaevi(3) ! Dipole polarizability tensor eigenvalues (imag) 1376 double precision ddotfile 1377 external ddotfile 1378 parameter (auXnm=45.563353d0) 1379 parameter (au2ang=5.29177249d-01) 1380 parameter (au2ang3=au2ang*au2ang*au2ang) 1381 logical nodezero 1382 character*4 irrepname 1383 character*3 axisname(3) ! Axis 1384 data axisname/'X','Y','Z'/ 1385 character*255 filename 1386 character*4 ir1filename(3) ! File name stub 1387 data ir1filename/'ir1x','ir1y','ir1z'/ 1388 character*4 ir2filename(3) ! File name stub 1389 data ir2filename/'ir2x','ir2y','ir1z'/ 1390 character*5 rr1filename(3) ! File name stub 1391 data rr1filename/'rr1x ','rr1y ','rr1z '/ 1392 character*5 rr2filename(3) ! File name stub 1393 data rr2filename/'rr2x ','rr2y ','rr2z '/ 1394c 1395 nodezero=(ga_nodeid().eq.0) 1396c 1397c CCSD-IR evaluation step 1398c 1399 cpu=-util_cpusec() 1400 wall=-util_wallsec() 1401 1402 omega = ifreq(omegacount) 1403 1404 do axisA = 1, 3 1405 do axisB = 1, axisA 1406 alpha(axisA,axisB)=0.0d0 1407 if (respaxis(axisA).and.respaxis(axisB)) then 1408 irrep_a=sym_abelian_axis(geom,axisA) 1409 irrep_b=sym_abelian_axis(geom,axisB) 1410 irrep_y=irrep_g 1411c 1412 call tce_filename('a0',filename) 1413 call createfile(filename,d_a0,1) ! size_a0 = 1 1414c 1415c if (nodezero) write(LuOut,*) "axisA = ",axisA 1416c if (nodezero) write(LuOut,*) "axisB = ",axisB 1417c 1418 alpha1=0.0d0 1419 alpha2=0.0d0 1420c 1421 irrep_d=irrep_a 1422 irrep_tr=irrep_b 1423c 1424c write(LuOut,*) "alpha_1 AB" 1425 call alpha_1(d_d1(axisA),d_a0,d_t1,d_t2, 1426 1 d_tr1(axisB+6),d_tr2(axisB+6), 1427 2 d_lambda1,d_lambda2, 1428 3 k_d1_offset(axisA),k_a0_offset, 1429 4 k_t1_offset,k_t2_offset, 1430 5 k_tr1_offset(axisB),k_tr2_offset(axisB), 1431 6 k_l1_offset,k_l2_offset) 1432c 1433 call reconcilefile(d_a0,1) 1434 call get_block(d_a0,alpha1,1,0) 1435 call tce_zero(d_a0,1) 1436c write(LuOut,*) "alpha_1 AB = ",alpha1 1437c 1438 if (axisA.eq.axisB) then 1439 alpha2=alpha1 1440 else 1441 irrep_d=irrep_b 1442 irrep_tr=irrep_a 1443c write(LuOut,*) "alpha_1 BA" 1444 call alpha_1(d_d1(axisB),d_a0,d_t1,d_t2, 1445 1 d_tr1(axisA+6),d_tr2(axisA+6), 1446 2 d_lambda1,d_lambda2, 1447 3 k_d1_offset(axisB),k_a0_offset, 1448 4 k_t1_offset,k_t2_offset, 1449 5 k_tr1_offset(axisA),k_tr2_offset(axisA), 1450 6 k_l1_offset,k_l2_offset) 1451c 1452 call reconcilefile(d_a0,1) 1453 call get_block(d_a0,alpha2,1,0) 1454 call tce_zero(d_a0,1) 1455 endif ! axisA.eq.axisB 1456c write(LuOut,*) "alpha_1 BA = ",alpha2 1457c 1458 alpha(axisA,axisB)=alpha(axisA,axisB)+alpha1+alpha2 1459c 1460 irrep_tra=irrep_a 1461 irrep_trb=irrep_b 1462c 1463 if (omega.ne.(0.0d0)) then 1464 alpha1=0.0d0 1465 alpha2=0.0d0 1466c write(LuOut,*) "alpha_2 A(+i) B(-i)" 1467 call alpha_2(d_f1,d_a0,d_t1,d_t2, 1468 1 d_tr1(axisA),d_tr2(axisA), 1469 2 d_tr1(axisB+3),d_tr2(axisB+3), 1470 3 d_v2,d_lambda1,d_lambda2,k_f1_offset,k_a0_offset, 1471 4 k_t1_offset,k_t2_offset, 1472 5 k_tr1_offset(axisA),k_tr2_offset(axisA), 1473 6 k_tr1_offset(axisB),k_tr2_offset(axisB),k_v2_offset, 1474 7 k_l1_offset,k_l2_offset) 1475 call reconcilefile(d_a0,1) 1476 call get_block(d_a0,alpha1,1,0) 1477 call tce_zero(d_a0,1) 1478c write(LuOut,*) "alpha_2 A(+i) B(-i) = ",alpha1 1479c 1480c write(LuOut,*) "alpha_2 A(-i) B(+i)" 1481 call alpha_2(d_f1,d_a0,d_t1,d_t2, 1482 1 d_tr1(axisA+3),d_tr2(axisA+3), 1483 2 d_tr1(axisB),d_tr2(axisB), 1484 3 d_v2,d_lambda1,d_lambda2,k_f1_offset,k_a0_offset, 1485 4 k_t1_offset,k_t2_offset, 1486 5 k_tr1_offset(axisA),k_tr2_offset(axisA), 1487 6 k_tr1_offset(axisB),k_tr2_offset(axisB),k_v2_offset, 1488 7 k_l1_offset,k_l2_offset) 1489 call reconcilefile(d_a0,1) 1490 call get_block(d_a0,alpha2,1,0) 1491 call tce_zero(d_a0,1) 1492c write(LuOut,*) "alpha_2 A(-i) B(+i) = ",alpha2 1493 elseif (omega.eq.(0.0d0)) then 1494 alpha1 = 0.0d0 1495 alpha2 = 0.0d0 1496 endif 1497c 1498 alpha3=0.0d0 1499c write(LuOut,*) "alpha_2 A(r) B(r)" 1500 call alpha_2(d_f1,d_a0,d_t1,d_t2, 1501 1 d_tr1(axisA+6),d_tr2(axisA+6), 1502 2 d_tr1(axisB+6),d_tr2(axisB+6), 1503 3 d_v2,d_lambda1,d_lambda2,k_f1_offset,k_a0_offset, 1504 4 k_t1_offset,k_t2_offset, 1505 5 k_tr1_offset(axisA),k_tr2_offset(axisA), 1506 6 k_tr1_offset(axisB),k_tr2_offset(axisB),k_v2_offset, 1507 7 k_l1_offset,k_l2_offset) 1508 call reconcilefile(d_a0,1) 1509 call get_block(d_a0,alpha3,1,0) 1510c write(LuOut,*) "alpha_2 A(r) B(r) = ",alpha3 1511 call deletefile(d_a0) 1512c 1513 alpha(axisA,axisB)=alpha(axisA,axisB)-alpha1-alpha2+alpha3 1514 alpha(axisA,axisB)=-1.0d0*alpha(axisA,axisB) 1515c 1516c write(LuOut,*) "alpha(axisA,axisB) = ",alpha(axisA,axisB) 1517c if (nodezero) write(LuOut,9020) cpu, wall 1518c 1519 endif ! respaxis(axis) 1520 enddo ! axisB loop 1521 enddo ! axisA loop 1522 cpu=cpu+util_cpusec() 1523 wall=wall+util_wallsec() 1524c 1525 do i = 1, 3 1526 do j = 1, i 1527 alphacopy(i,j)=alpha(i,j) 1528 alphacopy(j,i)=alpha(i,j) 1529 enddo 1530 enddo 1531c 1532 call hnd_diag(alphacopy,alphaevr,3,.false.,.false.) 1533c 1534 alphaiso = (alphaevr(1)+alphaevr(2)+alphaevr(3))/3.0d0 1535 alphaani = (alphaevr(1)-alphaevr(2))*(alphaevr(1)-alphaevr(2)) 1536 1 + (alphaevr(1)-alphaevr(3))*(alphaevr(1)-alphaevr(3)) 1537 2 + (alphaevr(2)-alphaevr(3))*(alphaevr(2)-alphaevr(3)) 1538 alphaani = dsqrt(0.5d0*alphaani) 1539c 1540 if ((.not.(respaxis(1).and.respaxis(2).and.respaxis(3))) 1541 1 .and.nodezero) write(LuOut,9911) 1542c 1543 if (nodezero) write(LuOut,9434) "CCSD Imaginary Response", 1544 1 ifreq(omegacount),auXnm/ifreq(omegacount), 1545 2 alpha(1,1),alpha(2,1),alpha(3,1), 1546 3 au2ang3*alpha(1,1),au2ang3*alpha(2,1),au2ang3*alpha(3,1), 1547 4 alpha(2,1),alpha(2,2),alpha(3,2), 1548 5 au2ang3*alpha(2,1),au2ang3*alpha(2,2),au2ang3*alpha(3,2), 1549 6 alpha(3,1),alpha(3,2),alpha(3,3), 1550 7 au2ang3*alpha(3,1),au2ang3*alpha(3,2),au2ang3*alpha(3,3), 1551 8 alphaevr(1),alphaevr(2),alphaevr(3), 1552 9 au2ang3*alphaevr(1),au2ang3*alphaevr(2),au2ang3*alphaevr(3), 1553 1 alphaiso,au2ang3*alphaiso,alphaani,au2ang3*alphaani 1554 if (nodezero) write(LuOut,9020) cpu, wall 1555 call util_flush(LuOut) 1556c 1557 ifreqval(omegacount,1) = alpha(1,1) 1558 ifreqval(omegacount,2) = alpha(2,2) 1559 ifreqval(omegacount,3) = alpha(3,3) 1560 ifreqval(omegacount,4) = alpha(2,1) 1561 ifreqval(omegacount,5) = alpha(3,1) 1562 ifreqval(omegacount,6) = alpha(3,2) 1563 ifreqval(omegacount,7) = alphaiso 1564 ifreqval(omegacount,8) = alphaani 1565c 1566 9020 format(1x,'Cpu & wall time / sec',2f15.1) 1567 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 1568 9100 format(1x,i4,2f18.13,2f8.1) 1569 9120 format(1x,A) 1570 9121 format(/,1x,A) 1571 9122 format(1x,A,i4) 1572 9400 format(/,1x,A,' iterations',/, 1573 1 1x,'---------------------------------------------',/ 1574 2 1x,'Iter Residuum Cpu Wall',/ 1575 3 1x,'---------------------------------------------') 1576 9410 format( 1577 1 1x,'---------------------------------------------',/ 1578 2 1x,'Iterations converged') 1579 9420 format(1x,i4,f25.13,2f8.1) 1580 9431 format(/,1x,'Frequency = ',f15.7,' / au') 1581 9434 format(/,1x,A,' polarizability / au ',/ 1582 1 1x,'Frequency = ',f15.7,' / au',/ 1583 1 1x,'Wavelength = ',f15.7,' / nm',/ 1584 3 1x,'-----------------------------------------------' 1585 3 ,'--------|-----------------------------------------------',/ 1586 2 1x,' atomic units (bohr^3) ' 1587 2 ,' | angstroms^3 ',/ 1588 2 1x,' X Y Z', 1589 2 1x,' | X Y Z',/ 1590 3 1x,'-----------------------------------------------' 1591 3 ,'--------|-----------------------------------------------',/ 1592 4 1x,'X ',3f15.7,3x,'|',3f15.7,/ 1593 5 1x,'Y ',3f15.7,3x,'|',3f15.7,/ 1594 6 1x,'Z ',3f15.7,3x,'|',3f15.7,/ 1595 3 1x,'-----------------------------------------------' 1596 3 ,'--------|-----------------------------------------------',/ 1597 6 1x,'Eigs = ',3f15.7,3x,'|',3f15.7,/ 1598 6 1x,'Isotropic = ',8x,1f15.7,3x,15x,'|',15x,1f15.7,/ 1599 6 1x,'Anisotropic = ',8x,1f15.7,3x,15x,'|',15x,1f15.7,/ 1600 3 1x,'-----------------------------------------------' 1601 3 ,'--------|-----------------------------------------------') 1602 9435 format(/,1x,A,' C6 coefficients ',/ 1603 1 1x,'--------------------------------',/ 1604 2 1x,'C6(XX) ',f15.7,/ 1605 3 1x,'C6(YY) ',f15.7,/ 1606 4 1x,'C6(ZZ) ',f15.7,/ 1607 5 1x,'C6(XY) ',f15.7,/ 1608 6 1x,'C6(XZ) ',f15.7,/ 1609 7 1x,'C6(YZ) ',f15.7,/ 1610 8 1x,'C6(AVG) ',f15.7,/ 1611 9 1x,'C6(ANI) ',f15.7,/ 1612 1 1x,'--------------------------------') 1613 9911 format(/,1x,'Warning: you have not solved ', 1614 1 'the response equations for all axes. ', 1615 2 'Please analyze the results carefully as ', 1616 3 'the average and anisotropic polarizabilities ', 1617 4 'are surely wrong.',/) 1618 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 1619 return 1620 end 1621 1622 1623 1624 1625 subroutine ccsd_ir_guess(d_d1,d_t1,d_t2,d_tr1,d_tr2, 1626 1 k_d1_offset,k_t1_offset,k_t2_offset, 1627 2 k_tr1_offset,k_tr2_offset, 1628 3 size_tr1,size_tr2,component,axis,omega) 1629 implicit none 1630#include "global.fh" 1631#include "mafdecls.fh" 1632#include "util.fh" 1633#include "errquit.fh" 1634#include "stdio.fh" 1635#include "tce.fh" 1636#include "tce_main.fh" 1637#include "tce_prop.fh" 1638c 1639 integer axis,component 1640 integer d_d1(3),d_t1,d_t2 1641 integer d_tr1(9),d_tr2(9),d_rr1(3),d_rr2(3) 1642 integer k_d1_offset(3),k_t1_offset,k_t2_offset 1643 integer k_tr1_offset(3),k_tr2_offset(3) 1644 integer size_tr1(3),size_tr2(3) 1645 double precision omega 1646 double precision cpu, wall 1647 logical nodezero 1648 character*4 irrepname 1649 character*3 axisname(3) ! Axis 1650 data axisname/'X','Y','Z'/ 1651 character*255 filename 1652 character*5 rr1filename(3) ! File name stub 1653 data rr1filename/'rr1x ','rr1y ','rr1z '/ 1654 character*5 rr2filename(3) ! File name stub 1655 data rr2filename/'rr2x ','rr2y ','rr2z '/ 1656c 1657 nodezero=(ga_nodeid().eq.0) 1658c 1659 if (component.eq.0) then 1660c 1661 if (nodezero) write(6,9121) 'Initial guess x = b/Adiag' 1662c 1663c 0. Create dummies 1664c 1665 call tce_filename(rr1filename(axis),filename) 1666 call createfile(filename,d_rr1(axis),size_tr1(axis)) 1667 call tce_zero(d_rr1(axis),size_tr1(axis)) 1668c 1669 call tce_filename(rr2filename(axis),filename) 1670 call createfile(filename,d_rr2(axis),size_tr2(axis)) 1671 call tce_zero(d_rr2(axis),size_tr2(axis)) 1672c 1673c 1. Form b 1674c 1675 call ccsd_o1(d_rr1(axis),d_d1(axis),d_t1,d_t2, 1676 1 k_tr1_offset(axis),k_d1_offset(axis), 1677 2 k_t1_offset,k_t2_offset) 1678 call ccsd_o2(d_rr2(axis),d_d1(axis),d_t1,d_t2, 1679 1 k_tr2_offset(axis),k_d1_offset(axis), 1680 2 k_t1_offset,k_t2_offset) 1681c 1682c 2. Hit with preconditioner 1683c 1684 call tce_jacobi_ir1(d_rr1(axis),d_tr1(axis+6), 1685 1 k_tr1_offset(axis),0.0d0,0.0d0,1.0d0) 1686 call tce_jacobi_ir2(d_rr2(axis),d_tr2(axis+6), 1687 1 k_tr2_offset(axis),0.0d0,0.0d0,1.0d0) 1688c 1689c 3. Delete dummies 1690c 1691 call deletefile(d_rr1(axis)) 1692 call deletefile(d_rr2(axis)) 1693c 1694 elseif (component.eq.1) then 1695c 1696c Initial guess x_I = omega*x_R/(Adiag-omega) 1697 if (nodezero) write(6,9121) 1698 1 'Initial guess x_I = omega*x_R/Adiag' 1699c 1700 if (omega.eq.(0.0d0)) then 1701 call tce_zero(d_tr1(axis+0),size_tr1(axis)) 1702 call tce_zero(d_tr2(axis+0),size_tr2(axis)) 1703 call tce_zero(d_tr1(axis+3),size_tr1(axis)) 1704 call tce_zero(d_tr2(axis+3),size_tr2(axis)) 1705 return 1706 endif 1707c 1708c 0. Create dummies 1709c 1710 call tce_filename(rr1filename(axis),filename) 1711 call createfile(filename,d_rr1(axis),size_tr1(axis)) 1712 call tce_zero(d_rr1(axis),size_tr1(axis)) 1713c 1714 call tce_filename(rr2filename(axis),filename) 1715 call createfile(filename,d_rr2(axis),size_tr2(axis)) 1716 call tce_zero(d_rr2(axis),size_tr2(axis)) 1717c 1718c 1. Form w_I = omega*x_R 1719c 1720 call daxfile(1,omega,d_tr1(axis+6), 1721 1 d_rr1(axis),size_tr1(axis)) 1722 call daxfile(1,omega,d_tr2(axis+6), 1723 1 d_rr2(axis),size_tr2(axis)) 1724c 1725c 3. Hit with preconditioner 1726c 1727 if (omega.gt.(0.0d0)) then 1728 call tce_jacobi_ir1(d_rr1(axis),d_tr1(axis+0), 1729 1 k_tr1_offset(axis),0.0d0,0.0d0,1.0d0) 1730 call tce_jacobi_ir2(d_rr2(axis),d_tr2(axis+0), 1731 1 k_tr2_offset(axis),0.0d0,0.0d0,1.0d0) 1732 elseif (omega.lt.(0.0d0)) then 1733 call tce_jacobi_ir1(d_rr1(axis),d_tr1(axis+3), 1734 1 k_tr1_offset(axis),0.0d0,0.0d0,1.0d0) 1735 call tce_jacobi_ir2(d_rr2(axis),d_tr2(axis+3), 1736 1 k_tr2_offset(axis),0.0d0,0.0d0,1.0d0) 1737 endif 1738c 1739c 4. Delete dummies 1740c 1741 call deletefile(d_rr1(axis)) 1742 call deletefile(d_rr2(axis)) 1743c 1744 endif ! component 1745c 1746 9020 format(1x,'Cpu & wall time / sec',2f15.1) 1747 9082 format(1x,'amplitude norm of ',A9,' = ',f25.15) 1748 9100 format(1x,i4,2f18.13,2f8.1) 1749 9120 format(1x,A) 1750 9121 format(/,1x,A) 1751 9122 format(1x,A,i4) 1752 9420 format(1x,i4,f25.13,2f8.1) 1753 9431 format(/,1x,'Frequency = ',f15.7,' / au') 1754 9440 format(1x,A3,' axis ( ',A4,'symmetry)') 1755 return 1756 end 1757 1758#ifdef DENOM_EXPER 1759c==================================================================== 1760c DENOMINATOR EXPERIMENT 1761c==================================================================== 1762 call tce_filename('denom1',filename) 1763 call createfile(filename,d_denom1,size_t1) 1764 call tce_zero(d_denom1,size_t1) 1765c 1766 call tce_filename('denom2',filename) 1767 call createfile(filename,d_denom2,size_t2) 1768 call tce_zero(d_denom2,size_t2) 1769c 1770 call tce_make_denom1(d_denom1,k_t1_offset,irrep_t, 1771 1 1,0.0d0,0.0d0) ! denom_power, omega, shift 1772 call tce_make_denom2(d_denom2,k_t2_offset,irrep_t, 1773 1 1,0.0d0,0.0d0) ! denom_power, omega, shift 1774c 1775 call tce_print_x1(d_denom1,k_t1_offset,1.0d-6,irrep_t) 1776 call tce_print_x2(d_denom2,k_t2_offset,1.0d-6,irrep_t) 1777c 1778 call deletefile(d_denom1) 1779 call deletefile(d_denom2) 1780c 1781 do axis = 1, 3 1782 do i = 1, 2 1783 if (respaxis(axis)) then 1784c 1785 if(nodezero) write(6,*) '=================================' 1786 if(nodezero) write(6,*) 'denom_power = ',i 1787c 1788 irrep_d=sym_abelian_axis(geom,axis) 1789 call sym_irrepname(geom,irrep_d+1,irrepname) 1790 if (nodezero.and.util_print('mod1',print_default)) then 1791 write(LuOut,*) 1792 write(LuOut,9440) axisname(axis),irrepname 1793 endif 1794c 1795 call tce_filename('denom1',filename) 1796 call createfile(filename,d_denom1,size_tr1(axis)) 1797 call tce_zero(d_denom1,size_tr1(axis)) 1798c 1799 call tce_filename('denom2',filename) 1800 call createfile(filename,d_denom2,size_tr2(axis)) 1801 call tce_zero(d_denom2,size_tr2(axis)) 1802c 1803 call tce_make_denom1(d_denom1,k_tr1_offset(axis),irrep_d, 1804 1 i,0.0d0,0.0d0) ! denom_power, omega, shift 1805 call tce_make_denom2(d_denom2,k_tr2_offset(axis),irrep_d, 1806 1 i,0.0d0,0.0d0) ! denom_power, omega, shift 1807c 1808 call tce_print_x1(d_denom1,k_tr1_offset(axis),1.0d-6, 1809 1 irrep_d) 1810 call tce_print_x2(d_denom2,k_tr2_offset(axis),1.0d-6, 1811 2 irrep_d) 1812c 1813 call deletefile(d_denom1) 1814 call deletefile(d_denom2) 1815c 1816 endif 1817 enddo 1818 enddo 1819c 1820 call errquit('tce_energy: MANUAL STOP',0,CALC_ERR) 1821c==================================================================== 1822#endif 1823c 1824c $Id$ 1825