1C $Id$ 2************************************************************************ 3c:tex-\subsection{rel\_onel} 4c:tex-This routine generates the modified one-electron integrals for 5c:tex-a relativistic basis set. These are the modified kinetic energy, 6c:tex-the modified potential energy and the modified overlap, 7c:tex-\begin{eqnarray} 8c:tex-&& \tilde{T}_{ab} = T_{ab}^{LS} + T_{ab}^{SL} - T_{ab}^{SS} 9c:tex- \nonumber \\ 10c:tex-&& \tilde{V}^{sf}_{ab} = V_{ab}^{LL} + {{\alpha^2}\over{4}} 11c:tex- \nabla_A\cdot\nabla_B V_{ab}^{SS} 12c:tex-&& \tilde{V}^{so}_{ab} = V_{ab}^{LL} + {{\alpha^2}\over{4}} 13c:tex- \nabla_A\times\nabla_B V_{ab}^{SS} 14c:tex- \nonumber \\ 15c:tex-&& \tilde{S}_{ab} = S_{ab}^{LL} + {{\alpha^2}\over{2}} T_{ab}^{SS} 16c:tex- \nonumber 17c:tex-\end{eqnarray} 18c:tex- 19c:tex-\noindent Author: K. G. Dyall 20c:tex- 21c:tex-{\it Syntax:} 22c:tex-\begin{verbatim} 23 subroutine rel_onel ( 24 & Axyz,zeta_A,coefL_A,coefS_A,n_prim_A,n_cont_A,l_A, 25 & Bxyz,zeta_B,coefL_B,coefS_B,n_prim_B,n_cont_B,l_B, 26 & Cxyz,zan,exinv,nat,S,T,V,lstv,doStil,doTtil,doVtil, 27 & canAB,do_nw,do_hnd,nonrel,DryRun,scr,lscr,ibug,ntyp) 28c:tex-\end{verbatim} 29 implicit none 30#include "stdio.fh" 31#include "rel_consts.fh" 32#include "errquit.fh" 33* 34c:tex-{\it Argument list:} 35c:tex-\begin{verbatim} 36 integer n_prim_A ! [input] num. prims in shell A 37 integer n_cont_A ! [input] num general conts in shell A 38 integer l_A ! [input] angular momentum of shell A 39 integer n_prim_B ! [input] num. prims in shell B 40 integer n_cont_B ! [input] num general conts in shell B 41 integer l_B ! [input] angular momentum of shell B 42 integer nat ! [input] number of atoms 43 integer lscr ! [input] size of scratch array 44 integer lstv ! [input] size of any integral buffer 45 integer ibug ! [input] debug variable 46 integer ntyp ! [input] potential energy integral type 47 double precision Axyz(3) ! [input] position of center A 48 double precision zeta_A(n_prim_A) ! [input] exponents of shell A 49 double precision coefL_A(n_prim_A,n_cont_A) ! [input] A large coeffs 50 double precision coefS_A(n_prim_A,n_cont_A) ! [input] A small coeffs 51 double precision Bxyz(3) ! [input] position of center B 52 double precision zeta_B(n_prim_B) ! [input] exponents of shell B 53 double precision coefL_B(n_prim_B,n_cont_B) ! [input] B large coeffs 54 double precision coefS_B(n_prim_B,n_cont_B) ! [input] B small coeffs 55 double precision Cxyz(3,nat) ! [input] all atom positions 56 double precision zan(nat) ! [input] charges on all atoms 57 double precision exinv(nat) ! [input] inverse nuclear exponents 58 double precision scr(lscr) ! [scratch] scratch buffers 59 double precision S(lstv) ! [output] overlap integrals 60 double precision T(lstv) ! [output] kinetic energy integrals 61 double precision V(lstv,ntyp) ! [output] potential integrals 62 logical doStil ! [input] compute modified overlap (True/False) 63 logical doTtil ! [input] compute modified kinetic (True/False) 64 logical doVtil ! [input] compute modified potential (True/False) 65 logical canAB ! [input] compute only canonical ints (false only) 66 logical do_nw ! [input] can do NW integrals 67 logical do_hnd ! [input] can do HONDO integrals 68 logical nonrel ! [input] true if either centre is nonrelativistic 69 logical DryRun ! [input] true means only compute required memory 70c:tex-\end{verbatim} 71c:tex-See rel_pot for a description of the allowed values of ibug and ntyp 72c:tex- 73c:tex-{\it Subroutines called:} int\_hf1sp, rel\_pot, daxpy 74* 75 integer n_cart_a ! cartesian components of shell A 76 integer n_cart_b ! cartesian components of shell B 77 integer n_cart_ab ! n_cart_a*n_cart_b 78 integer n_cont_ab ! n_cont_a*n_cont_b 79 integer n_all_b ! n_cart_b*n_cont_b 80 integer n_all_a ! n_cart_a*n_cont_a 81 integer n_ab ! number of integrals 82 integer n_cartp_a ! cartesian components for l_A+1 83 integer n_cartp_b ! cartesian components for l_B+1 84 integer n_allp_b ! n_cartp_b*n_cont_b 85 integer n_allp_a ! n_cartp_a*n_cont_a 86 integer n_cartm_a ! cartesian components for l_A-1 87 integer n_cartm_b ! cartesian components for l_B-1 88 integer n_allm_b ! n_cartm_b*n_cont_b 89 integer n_allm_a ! n_cartm_a*n_cont_a 90 integer n_intpp ! number of integrals for l_A+1,l_B+1 91 integer n_intpm ! number of integrals for l_A-1,l_B+1 92 integer n_intmp ! number of integrals for l_A+1,l_B-1 93 integer n_intmm ! number of integrals for l_A-1,l_B-1 94 integer i_xca ! address in scr of exp*coef for shell A 95 integer i_xcb ! address in scr of exp*coef for shell B 96 integer i_pp ! address in scr of integrals for l_A+1,l_B+1 97 integer i_pm ! address in scr of integrals for l_A-1,l_B+1 98 integer i_mp ! address in scr of integrals for l_A+1,l_B-1 99 integer i_mm ! address in scr of integrals for l_A-1,l_B-1 100 integer i_scr ! address of free space in scr 101 integer memscr ! free space in scr 102 integer max_mem ! maximum memory used 103 integer i,j,k ! loop indices etc. 104 double precision one ! Obvious! 105 parameter (one = 1.0D0) 106* 107 logical debug_gen ! do general debug printing 108 logical debug_addresses ! do address debug printing 109 logical debug_arrays ! do array debug printing 110 logical doS ! compute overlap (True/False) 111 logical doT ! compute kinetic (True/False) 112 logical doV ! compute potential (True/False) 113 character*12 pot_type(4) 114 data pot_type 115 & /' Scalar','z spin-orbit','y spin-orbit','x spin-orbit'/ 116* 117 debug_gen = ibug .gt. 0 118 debug_addresses = mod(ibug,2) .eq. 1 119 debug_arrays = (mod(ibug,10)/2 .eq. 1) .and. .not.DryRun 120 max_mem = 0 121* 122 if (.not.(do_nw .or. do_hnd)) call errquit 123 & ('rel_onel: can''t do NW or HONDO integrals',99, INT_ERR) 124* 125 if (debug_gen) then 126 write (LuOut,*) 127 write (LuOut,*) 'rel_onel: ibug =',ibug 128 write (LuOut,*) 'l_A,n_prim_A,n_cont_A',l_A,n_prim_A,n_cont_A 129 write (LuOut,*) 'l_B,n_prim_B,n_cont_B',l_B,n_prim_B,n_cont_B 130 if (do_hnd) then 131 write (LuOut,*) 'Using HONDO integrals' 132 else if (do_nw) then 133 write (LuOut,*) 'Using NWChem integrals' 134 end if 135 end if 136* 137 n_cart_a = (l_a+1)*(l_a+2)/2 138 n_cart_b = (l_b+1)*(l_b+2)/2 139 n_cart_ab = n_cart_a*n_cart_b 140 n_cont_ab = n_cont_a*n_cont_b 141 n_all_a = n_cart_a*n_cont_a 142 n_all_b = n_cart_b*n_cont_b 143 n_ab = n_cart_ab*n_cont_ab 144 if (lstv .lt. n_ab .and. .not.DryRun) call errquit ( 145 & 'Integral buffer length too small in rel_onel',99, 146 & MEM_ERR) 147 if (debug_addresses) then 148 write (LuOut,*) 'n_cart_a',n_cart_a 149 write (LuOut,*) 'n_cart_b',n_cart_b 150 write (LuOut,*) 'n_cart_ab',n_cart_ab 151 write (LuOut,*) 'n_cont_ab',n_cont_ab 152 write (LuOut,*) 'n_all_a',n_all_a 153 write (LuOut,*) 'n_all_b',n_all_b 154 write (LuOut,*) 'n_ab',n_ab 155 end if 156 if (debug_arrays) then 157 call ecp_matpr (coefL_A,1,n_prim_a,1,n_cont_a, 158 & 1,n_prim_a,1,n_cont_a,'L coef A','E',120,6) 159 call ecp_matpr (coefS_A,1,n_prim_a,1,n_cont_a, 160 & 1,n_prim_a,1,n_cont_a,'S coef A','E',120,6) 161 call ecp_matpr (coefL_B,1,n_prim_b,1,n_cont_b, 162 & 1,n_prim_b,1,n_cont_b,'L coef B','E',120,6) 163 call ecp_matpr (coefS_B,1,n_prim_b,1,n_cont_b, 164 & 1,n_prim_b,1,n_cont_b,'S coef B','E',120,6) 165 end if 166* 167* Calculate large component overlap and nuclear attraction integrals, 168* and kinetic energy integrals if center is nonrelativistic. 169* 170 doS = doStil 171 doT = doTtil.and.nonrel 172 doV = doVtil 173 memscr = lscr 174 if (do_hnd) then 175 call hnd_stvint( 176 & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A, 177 & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B, 178 & Cxyz,zan,nat,S,T,V,n_ab,doS,doT,doV,scr,memscr) 179 else if (do_nw) then 180 call hf1( 181 & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A, 182 & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B, 183 & Cxyz,zan,exinv,nat,S,T,V,n_ab,doS,doT,doV,canAB, 184 & DryRun,scr,memscr) 185 end if 186 if (debug_arrays) then 187 if (doS) call ecp_matpr(S,1,n_all_b,1,n_all_a, 188 & 1,n_all_b,1,n_all_a,'LL overlap','E',120,6) 189 if (doT) call ecp_matpr(T,1,n_all_b,1,n_all_a, 190 & 1,n_all_b,1,n_all_a,'LL kinetic','E',120,6) 191 if (doV) call ecp_matpr(V,1,n_all_b,1,n_all_a, 192 & 1,n_all_b,1,n_all_a,'LL potential','E',120,6) 193 end if 194 if (DryRun) max_mem = max(max_mem,memscr) 195 if (nonrel) return 196* 197* Calculate kinetic energy integrals, correction to overlaps 198* 199 if (doTtil) then 200 doS = .false. 201 doT = .true. 202 doV = .false. 203 memscr = lscr-n_ab 204 if (do_hnd) then 205 call hnd_stvint( 206 & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A, 207 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B, 208 & Cxyz,zan,nat,S,T,V,n_ab,doS,doT,doV,scr(n_ab+1),memscr) 209 else if (do_nw) then 210 call hf1( 211 & Axyz,zeta_A,coefL_A,n_prim_A,n_cont_A,l_A, 212 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B, 213 & Cxyz,zan,exinv,nat,S,T,V,n_ab,doS,doT,doV,canAB, 214 & DryRun,scr(n_ab+1),memscr) 215 end if 216 if (DryRun) then 217 max_mem = max(max_mem,memscr) 218 else if (debug_arrays) then 219 call ecp_matpr (T,1,n_all_b,1,n_all_a, 220 & 1,n_all_b,1,n_all_a,'LS kinetic','E',120,6) 221 end if 222 memscr = lscr-n_ab 223 if (do_hnd) then 224 call hnd_stvint( 225 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A, 226 & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B, 227 & Cxyz,zan,nat,S,scr,V,n_ab,doS,doT,doV,scr(n_ab+1),memscr) 228 else if (do_nw) then 229 call hf1( 230 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A, 231 & Bxyz,zeta_B,coefL_B,n_prim_B,n_cont_B,l_B, 232 & Cxyz,zan,exinv,nat,S,scr,V,n_ab,doS,doT,doV,canAB, 233 & DryRun,scr(n_ab+1),memscr) 234 end if 235 if (DryRun) then 236 max_mem = max(max_mem,memscr+n_ab) 237 else 238 if (debug_arrays) call ecp_matpr (scr,1,n_all_b,1,n_all_a, 239 & 1,n_all_b,1,n_all_a,'SL kinetic','E',120,6) 240 call daxpy (n_ab,one,scr,1,T,1) 241 if (debug_arrays) call ecp_matpr (T,1,n_all_b,1,n_all_a, 242 & 1,n_all_b,1,n_all_a,'LS+SL kinetic','E',120,6) 243 end if 244 end if 245 if (doStil .or. doTtil) then 246 doS = .false. 247 doT = .true. 248 doV = .false. 249 memscr = lscr-n_ab 250 if (do_hnd) then 251 call hnd_stvint( 252 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A, 253 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B, 254 & Cxyz,zan,nat,S,scr,V,n_ab,doS,doT,doV,scr(n_ab+1),memscr) 255 else if (do_nw) then 256 call hf1( 257 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A, 258 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B, 259 & Cxyz,zan,exinv,nat,S,scr,V,n_ab,doS,doT,doV,canAB, 260 & DryRun,scr(n_ab+1),memscr) 261 end if 262 if (DryRun) then 263 max_mem = max(max_mem,memscr+n_ab) 264 else 265 if (debug_arrays) call ecp_matpr (scr,1,n_all_b,1,n_all_a, 266 & 1,n_all_b,1,n_all_a,'SS kinetic','E',120,6) 267 if (doTtil) then 268 call daxpy (n_ab,-one,scr,1,T,1) 269 if (debug_arrays) call ecp_matpr (T,1,n_all_b,1,n_all_a, 270 & 1,n_all_b,1,n_all_a,'LS+SL-SS kinetic','E',120,6) 271 end if 272 if (doStil) then 273 call daxpy (n_ab,halsq,scr,1,S,1) 274 if (debug_arrays) call ecp_matpr (S,1,n_all_b,1,n_all_a, 275 & 1,n_all_b,1,n_all_a,'LL+SS overlap','E',120,6) 276 end if 277 end if 278 end if 279* 280* Generate small component potential arrays 281* 282* 283* Set up pointers to scratch space for coefficients multiplied by 284* exponents and for integrals with shifted l values 285* 286 if (.not.doVtil) return 287* 288 n_cartp_a = n_cart_a+l_A+2 289 n_cartp_b = n_cart_b+l_B+2 290 n_cartm_a = n_cart_a-l_A-1 291 n_cartm_b = n_cart_b-l_B-1 292 n_allp_a = n_cartp_a*n_cont_A 293 n_allp_b = n_cartp_b*n_cont_B 294 n_allm_a = n_cartm_a*n_cont_A 295 n_allm_b = n_cartm_b*n_cont_B 296 n_intpp = n_allp_a*n_allp_b 297 n_intpm = n_allm_a*n_allp_b 298 n_intmp = n_allp_a*n_allm_b 299 n_intmm = n_allm_a*n_allm_b 300 i_xca = 1 301 i_xcb = i_xca+n_prim_A*n_cont_A 302 i_pp = i_xcb+n_prim_B*n_cont_B 303 i_pm = i_pp+n_intpp 304 i_mp = i_pm+n_intpm 305 i_mm = i_mp+n_intmp 306 i_scr = max(i_xca+n_ab*ntyp,i_mm+n_intmm) 307* 308 if (debug_addresses) then 309 write (LuOut,*) 'n_cartp_a',n_cartp_a 310 write (LuOut,*) 'n_cartp_b',n_cartp_b 311 write (LuOut,*) 'n_cartm_a',n_cartm_a 312 write (LuOut,*) 'n_cartm_b',n_cartm_b 313 write (LuOut,*) 'n_intpp',n_intpp 314 write (LuOut,*) 'n_intpm',n_intpm 315 write (LuOut,*) 'n_intmp',n_intmp 316 write (LuOut,*) 'n_intmm',n_intmm 317 write (LuOut,*) 'i_xca,i_xcb',i_xca,i_xcb 318 write (LuOut,*) 'i_pp,i_pm,i_mp,i_mm',i_pp,i_pm,i_mp,i_mm 319 write (LuOut,*) 'i_scr',i_scr 320 end if 321* 322* Set up coefficients multiplied by exponents 323* 324 memscr = lscr-i_scr+1 325 if (.not.DryRun) then 326 if (memscr .lt. 0) call errquit ( 327 & 'Insufficient scratch memory in rel_onel',99, 328 & MEM_ERR) 329 k = i_xca-1 330 do j = 1,n_cont_A 331 do i = 1,n_prim_A 332 scr(k+i) = zeta_A(i)*coefS_A(i,j) 333 end do 334 k = k+n_prim_A 335 end do 336 k = i_xcb-1 337 do j = 1,n_cont_B 338 do i = 1,n_prim_B 339 scr(k+i) = zeta_B(i)*coefS_B(i,j) 340 end do 341 k = k+n_prim_B 342 end do 343 end if 344 doS = .false. 345 doT = .false. 346 doV = .true. 347* 348* Calculate integrals for l_A+1, l_B+1 349* 350 if (do_hnd) then 351 call hnd_stvint( 352 & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1, 353 & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1, 354 & Cxyz,zan,nat,S,T,scr(i_pp),n_intpp,doS,doT,doV, 355 & scr(i_scr),memscr) 356 else if (do_nw) then 357 call hf1( 358 & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1, 359 & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1, 360 & Cxyz,zan,exinv,nat,S,T,scr(i_pp),n_intpp,doS,doT,doV, 361 & canAB,DryRun,scr(i_scr),memscr) 362 end if 363 if (DryRun) then 364 max_mem = max(max_mem,i_scr+memscr-1) 365 if (debug_addresses) 366 & write (luout,*) '++:memscr,max_mem = ',memscr,max_mem 367 memscr = lscr-i_scr+1 368 else if (debug_arrays) then 369 call ecp_matpr (scr(i_pp),1,n_allp_b,1,n_allp_a, 370 & 1,n_allp_b,1,n_allp_a,'V(la+1,lb+1)','E',120,6) 371 end if 372* 373* Calculate integrals for l_A-1, l_B+1 374* 375 if (l_A .gt. 0) then 376 if (do_hnd) then 377 call hnd_stvint( 378 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1, 379 & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1, 380 & Cxyz,zan,nat,S,T,scr(i_pm),n_intpm,doS,doT,doV, 381 & scr(i_scr),memscr) 382 else if (do_nw) then 383 call hf1( 384 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1, 385 & Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1, 386 & Cxyz,zan,exinv,nat,S,T,scr(i_pm),n_intpm,doS,doT,doV, 387 & canAB,DryRun,scr(i_scr),memscr) 388 end if 389 if (DryRun) then 390 max_mem = max(max_mem,i_scr+memscr-1) 391 if (debug_addresses) 392 & write (luout,*) '-+:memscr,max_mem = ',memscr,max_mem 393 memscr = lscr-i_scr+1 394 else if (debug_arrays) then 395 call ecp_matpr (scr(i_pm),1,n_allp_b,1,n_allm_a, 396 & 1,n_allp_b,1,n_allp_a,'V(la-1,lb+1)','E',120,6) 397 end if 398 end if 399* 400* Calculate integrals for l_A+1, l_B-1 401* 402 if (l_B .gt. 0) then 403 if (do_hnd) then 404 call hnd_stvint( 405 & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1, 406 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1, 407 & Cxyz,zan,nat,S,T,scr(i_mp),n_intmp,doS,doT,doV, 408 & scr(i_scr),memscr) 409 else if (do_nw) then 410 call hf1( 411 & Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1, 412 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1, 413 & Cxyz,zan,exinv,nat,S,T,scr(i_mp),n_intmp,doS,doT,doV, 414 & canAB,DryRun,scr(i_scr),memscr) 415 end if 416 if (DryRun) then 417 max_mem = max(max_mem,i_scr+memscr-1) 418 if (debug_addresses) 419 & write (luout,*) '+-:memscr,max_mem = ',memscr,max_mem 420 memscr = lscr-i_scr+1 421 else if (debug_arrays) then 422 call ecp_matpr (scr(i_mp),1,n_allm_b,1,n_allp_a, 423 & 1,n_allm_b,1,n_allp_a,'V(la+1,lb-1)','E',120,6) 424 end if 425* 426* Calculate integrals for l_A-1, l_B-1 427* 428 if (l_A .gt. 0) then 429 if (do_hnd) then 430 call hnd_stvint( 431 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1, 432 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1, 433 & Cxyz,zan,nat,S,T,scr(i_mm),n_intmm,doS,doT,doV, 434 & scr(i_scr),memscr) 435 else if (do_nw) then 436 call hf1( 437 & Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1, 438 & Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1, 439 & Cxyz,zan,exinv,nat,S,T,scr(i_mm),n_intmm,doS,doT,doV, 440 & canAB,DryRun,scr(i_scr),memscr) 441 end if 442 if (DryRun) then 443 max_mem = max(max_mem,i_scr+memscr-1) 444 if (debug_addresses) 445 & write (luout,*) '--:memscr,max_mem = ',memscr,max_mem 446 memscr = lscr-i_scr+1 447 else if (debug_arrays) then 448 call ecp_matpr (scr(i_mm),1,n_allm_b,1,n_allm_a, 449 & 1,n_allm_b,1,n_allm_a,'V(la+1,lb-1)','E',120,6) 450 end if 451 end if 452 end if 453* 454* Compute the relativistic potential energy integrals 455* 456 call rel_pot (scr(i_pp),scr(i_pm),scr(i_mp),scr(i_mm), 457 & scr,n_ab,ntyp, 458 & l_A,n_cartp_a,n_cart_a,n_cartm_a,n_cont_A, 459 & l_B,n_cartp_b,n_cart_b,n_cartm_b,n_cont_B, 460 & DryRun,scr(i_scr),memscr,ibug/10) 461 if (DryRun) then 462 if (debug_addresses) write (luout,*) 'rel_pot:max_mem = ', 463 & max_mem 464 max_mem = max(max_mem,i_scr+memscr-1) 465 lscr = max_mem 466 if (debug_addresses) write (luout,*) 'max_mem = ',max_mem 467 else 468 i = 1 469 do j = 1,ntyp 470 if (debug_arrays) then 471 write (LuOut,'(//2A)') pot_type(j),' potential' 472 call ecp_matpr (scr(i),1,n_all_b,1,n_all_a,1,n_all_b, 473 & 1,n_all_a,'SS potential','E',120,6) 474 end if 475 call daxpy (n_ab,qalsq,scr(i),1,V(1,j),1) 476 i = i+n_ab 477 if (debug_arrays) call ecp_matpr (V(1,j),1,n_all_b,1,n_all_a, 478 & 1,n_all_b,1,n_all_a,'Relativistic potential','E',120,6) 479 end do 480 end if 481* 482 return 483 end 484