1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 subroutine cspa(ioptlwf,iopt,natoms,no_u,no_l,lasto,isa, 9 . qa,rcoor,rh,cell,xa,nhmax,numh,listh,listhptr, 10 . maxnc,ncmax,nctmax,nfmax,nftmax,nhijmax,nbands, 11 . no_cl,nspin,Node) 12C ****************************************************************************** 13C This subroutine builds the Localized Wave Functions, centered 14C on ATOMS (searching for atoms within a cutoff radius rcoor), and 15C assigns a RANDOM initial guess to start the CG minimization. 16C 17C Criterion to build LWF's: 18C 1) Method of Kim et al: use more localized orbitals than 19C occupied orbitals. 20C We assign the minimum number of orbitals so that there 21C is place for more electrons than those in the system; 22C for instance: 23C H: 1 LWF 24C C,Si: 3 LWF's 25C N: 3 LWF's 26C O: 4 LWF's 27C ... 28C 2) Method of Ordejon et al: number of localized orbitals 29C equal to number of occupied orbitals. 30C For the initial assignment of LWF centers to atoms, atoms 31C with even number of electrons, n, get n/2 LWFs. Odd atoms 32C get (n+1)/2 and (n-1)/2 in an alternating sequence, ir order 33C of appearance (controlled by the input in the atomic coor block). 34C 35C Written by P.Ordejon, 1993. 36C Re-written by P.Ordejon, November'96. 37C Corrected by P.Ordejon, April'97, May'97 38C lmax, lmaxs and nzls erased from the input by DSP, Aug 1998. 39C Alternating sequence for odd species in Ordejon, by E.Artacho, Aug 2008. 40C ******************************* INPUT *************************************** 41C integer ioptlwf : Build LWF's according to: 42C 0 = Read blindly from disk 43C 1 = Functional of Kim et al. 44C 2 = Functional of Ordejon-Mauri 45C integer iopt : 0 = Find structure of sparse C matrix and 46C build initial guess 47C 1 = Just find structure of C matrix 48C integer natoms : Number of atoms 49C integer no_u : Number of basis orbitals 50C integer lasto(0:natoms) : Index of last orbital of each atom 51C integer isa(natoms) : Species index of each atom 52C real*8 qa(natoms) : Neutral atom charge 53C real*8 rcoor : Cutoff radius of Localized Wave Functions 54C real*8 rh : Maximum cutoff radius of Hamiltonian matrix 55C real*8 cell(3,3) : Supercell vectors 56C real*8 xa(3,natoms) : Atomic coordinates 57C integer maxnc : First dimension of C matrix, and maximum 58C number of nonzero elements of each row of C 59C integer nspin : Number of spins 60C ****************************** OUTPUT ************************************** 61C real*8 c(ncmax,no_u) : Initial guess for sparse coefficients 62C of LWF's (only if iopt = 0) 63C integer numc(no_u) : Control vector of C matrix 64C integer listc(ncmax,no_u): Control vector of C matrix 65C integer ncmax : True value for ncmax, 66C If ncmax it is too small, then 67C c, numc and listc are NOT initialized!! 68C integer nctmax : Maximum number of nonzero elements 69C of eaxh column of C 70C integer nfmax : Maximum number of nonzero elements 71C of each row of F 72C integer nftmax : Maximum number of nonzero elements 73C of eaxh column of F 74C integer nhijmax : Maximum number of nonzero elements 75C of each row of Hij 76C integer nbands : Number of LWF's 77C **************************************************************************** 78 79 use precision, only: dp 80 use alloc, only: re_alloc, de_alloc 81 use on_main, only: numc, listc, c, cold, listcold 82 use on_main, only: ncg2l, ncl2g, nct2p, ncp2t 83 use neighbour, only: jan, r2ij, xij, mneighb, maxnna, 84 & reset_neighbour_arrays 85 use sys, only : die 86 use spatial, only : nL2G 87 88 implicit none 89 90 integer, intent(in) :: 91 . iopt,ioptlwf,natoms,no_u,nhmax,Node,no_l,nspin 92 93 integer, intent(inout) :: maxnc 94 95 integer, intent(out) :: 96 $ nbands,ncmax,nctmax, nfmax,nftmax,nhijmax, 97 $ no_cl 98 99 integer, intent(in) :: 100 . isa(natoms),lasto(0:natoms), 101 . numh(no_l),listh(nhmax),listhptr(no_l) 102 103 real(dp), intent(in) :: 104 . cell(3,3),qa(natoms),rcoor,rh,xa(3,natoms) 105 106C Internal variables ....................................................... 107 108 integer, dimension(:), pointer :: alist, ilist, numft 109 integer, pointer :: indexloc(:) => null() 110 111 integer 112 . i,ia,imu,in,index,indexa,indexb,indexi,indexj,iorb,is,j,ja, 113 . jj,k,mu,nct,nelectr,nf,nm,norb,nqtot,nu,numloc,nhij,indmu, 114 . mul, mull, ind, nelectr1, nelectr2 115 116 integer, save :: iseed = 17 117 integer, save :: nna = 200 118 119 logical, dimension(:), pointer :: lneeded 120 121 real(dp) :: 122 . cg,fact(2),qtot,r,rmax,rr(3),rrmod, 123 . snor,tiny,randomg,cgval 124 125 external :: randomg 126 127 logical, save :: firstcall = .true., secondodd = .false. 128 129 if (firstcall) then 130 131C Initialise random number generator 132 cgval = randomg(-iseed) 133 134 firstcall = .false. 135 136 endif 137 138 tiny = 1.d-10 139 140C Check that iopt is correct ................................................ 141 if (iopt .ne. 0 .and. iopt .ne. 1) then 142 call die('cspa: Wrong iopt in cspa') 143 endif 144 145 146C Allocate local arrays 147 nullify( alist, ilist, numft, lneeded ) 148 call re_alloc( alist, 1, natoms, 'alist', 'cspa' ) 149 call re_alloc( ilist, 1, no_u, name='ilist', routine='cspa' ) 150 call re_alloc( numft, 1, no_u, name='numft', routine='cspa' ) 151 call re_alloc( lneeded, 1, no_u, name='lneeded', 152 & routine='cspa' ) 153 154C Work out which basis functions are required in local C copy 155! These will be those which interact with the local orbitals 156 157 lneeded(1:no_u) = .false. 158 do mul = 1,no_l 159 ind = listhptr(mul) 160 lneeded(nL2G(mul,Node+1)) = .true. 161 do i = 1,numh(mul) 162 mu = listh(ind+i) 163 lneeded(mu) = .true. 164 enddo 165 enddo 166! 167! Now build up the indexes ncG2L and ncL2G 168! 169 no_cl = 0 170 ncG2L(1:no_u) = 0 171 ncT2P(1:no_u) = 0 172 ! First set of orbitals: those handled by 173 ! the local node 174 do mu = 1,no_l 175 no_cl = no_cl + 1 ! = mu 176 ncL2G(no_cl) = nL2G(mu,Node+1) 177 ncG2L(nL2G(mu,Node+1)) = no_cl ! = mu 178 enddo 179 ! Second set of orbitals: those interacting 180 ! with the orbitals handled by the local node 181 ! (and not already counted) 182 do mu = 1,no_u 183 if (lneeded(mu).and.ncG2L(mu).eq.0) then 184 no_cl = no_cl + 1 185 ncL2G(no_cl) = mu 186 ncG2L(mu) = no_cl 187 endif 188 enddo 189 do mul = 1,no_l 190 ! These are just identity mappings, 191 ! since the first no_l entries are 192 ! the same in the no_l and the no_cl sets 193 194 ! nct2p is just a filter: if >0, its argument 195 ! belongs to the (partial) set of indexes 196 ! that go from 1 to no_l. It could really be 197 ! a logical mask. It is always used as 198 ! if (nct2p(i)>0) then... 199 200 ! ncp2t is just the identity over 1:no_l, 201 ! where it is used in the rest of the program 202 203 ! Both arrays are dimensioned to no_u, even though 204 ! only the first no_l entries are used for ncP2T. 205 206 ncP2T(mul) = ncG2L(nL2G(mul,Node+1)) ! ncP2T(mul) = mul 207 ! Note also that the range of ncG2L is 1:no_cl, 208 ncT2P(ncG2L(nL2G(mul,Node+1))) = mul ! ncT2P(mul) = mul 209 enddo 210 211C Resize arrays that depend on no_cl 212 call re_alloc(listc,1,maxnc,1,no_cl, 213 . shrink=.false.,name='listc') 214 call re_alloc(listcold,1,maxnc,1,no_l, 215 . shrink=.false.,name='listcold') 216 call re_alloc(c,1,maxnc,1,no_cl,1,nspin, 217 . shrink=.false.,name='c') 218 call re_alloc(cold,1,maxnc,1,no_l,1,nspin, 219 . shrink=.false.,name='cold') 220 221C Initialize some stuff 222 223 do ia = 1,natoms 224 alist(ia) = 0 225 enddo 226 227 do mul = 1,no_cl 228 numc(mul) = 0 229 enddo 230 231 do mu = 1,no_u 232 numft(mu) = 0 233 enddo 234 235 if (iopt .eq. 0) then 236 do is = 1,nspin 237 do mu = 1,no_cl 238 do i = 1,maxnc 239 c(i,mu,is) = 0.0d0 240 enddo 241 enddo 242 enddo 243 endif 244 245 ncmax = 0 246 nctmax = 0 247 nfmax = 0 248 nftmax = 0 249 nhijmax = 0 250C ........................ 251 252C Calculate maximum length in unit cell ...................................... 253C determine maximum cell length 254 !AG: Possible bug: cell vectors are given by the columns of cell! 255 rmax = 0.0d0 256 do i = -1,1 257 do j = -1,1 258 do k = -1,1 259 !AG It should be 260 ! rr(1) = i*cell(1,1) + j*cell(1,2) + k*cell(1,3) 261 rr(1) = i*cell(1,1) + j*cell(2,1) + k*cell(3,1) 262 rr(2) = i*cell(1,2) + j*cell(2,2) + k*cell(3,2) 263 rr(3) = i*cell(1,3) + j*cell(2,3) + k*cell(3,3) 264 rrmod = sqrt( rr(1)**2 + rr(2)**2 + rr(3)**2 ) 265 if (rrmod .gt. rmax) rmax = rrmod 266 enddo 267 enddo 268 enddo 269C ........................ 270 271C Check that there is an integer number of electrons 272 qtot = 0.0d0 273 do ia = 1,natoms 274 qtot = qtot + qa(ia) 275 enddo 276 qtot = qtot + tiny 277 nqtot = nint(qtot) 278 if (abs(nqtot-qtot+tiny) .gt. 1e-3) then 279 if (Node.eq.0) then 280 write(6,*) 'cspa: Wrong total charge; non integer:',qtot 281 endif 282 call die() 283 endif 284C .................. 285 286C Build control vectors of sparse LWF 287C loop over the localized wave funcions (centered at atoms).................. 288 289C Initialize routine for neighbour search 290 ! Note that this will not update nna!! 291 call mneighb(cell,rcoor,natoms,xa,0,0,nna) 292 293C Allocate local memory 294 call re_alloc(indexloc,1,max(maxnna,natoms),name='indexloc') 295 296 index = 0 ! Counter for total number of LWF 297 loop_ia: do ia = 1,natoms 298 299 if (2.0*rcoor .lt. rmax) then 300 ! Localization size smaller than max cell size 301 ! Look for neighbours of atom ia 302 call mneighb(cell,rcoor,natoms,xa,ia,0,nna) 303 ! Reallocate in case nna has increased (Tight form) 304 call re_alloc( indexloc, 1, nna, 'indexloc', 'cspa' ) 305 else 306 ! Localization size greater than max cell size 307 ! We can treat all atoms in the cell as neighbors 308 ! Make sure we have enough space 309 ! (Cleaner alternative to initial allocation with 310 ! max(maxnna,natoms) 311 !call sizeup_neighbour_arrays(natoms) 312 !call re_alloc( indexloc, 1, natoms, name='indexloc') 313 nna = natoms 314 do jj = 1,natoms 315 jan(jj) = jj 316 enddo 317 endif 318 319 ! All procs have information about the number of LWFs 320 ! but they keep only the coefficients associated to 321 ! the orbitals they manage and to those that interact with 322 ! them 323 324 ! loop over LWF's centered on atom ia 325 call get_number_of_lwfs_on_atom() ! gets indexi, nelectr 326 loop_lwfs_on_ia: do indexb = 1,indexi 327 index = index + 1 ! global counter for number of LWFs 328 329c clear list of atoms considered within loc. range 330 do indexa = 1,nna 331 indexloc(indexa) = 0 332 enddo 333 numloc = 0 334 335C initialize stuff... 336 nct = 0 ! total number of coeffs for this LWF in this node 337 snor = 0.0d0 338 nf = 0 339 do nu = 1,no_u 340 ilist(nu) = 0 341 enddo 342 343C loop over the neighbors of ia within rcoor 344 loop_neighbor_atoms: do jj = 1,nna 345 ja = jan(jj) 346 347 ! Check if ja has already been included in current lwf 348 ! (an image atom, maybe?) 349 do indexa = 1,numloc 350 if (ja .eq. indexloc(indexa)) cycle loop_neighbor_atoms 351 enddo 352 numloc = numloc + 1 353 indexloc(numloc) = ja 354 355 ! Loop over orbitals of ja 356 norb = lasto(ja) - lasto(ja-1) 357 loop_orbs_on_neighbor: do iorb = 1,norb 358 mu = iorb + lasto(ja-1) 359 360 ! Get random number here for reproducibility over numbers of processors 361 get_random: do 362 cgval = (randomg(iseed) - 0.5d0) * 2.0d0 363 if (abs(cgval) .ge. 1d-5) exit get_random 364 enddo get_random 365 366 ! If this orbital (global index mu) is part of our 367 ! c-list, include it in the indexes 368 369 !! alternative: 370 !! if (ncG2L(mu) == 0) cycle loop_orbs_on_neighbor 371 372 if (ncG2L(mu).ne.0) then 373 mul = ncG2L(mu) 374 nm = numc(mul) 375 numc(mul) = nm + 1 376 if (numc(mul) .gt. maxnc) then 377 ! Reallocate all arrays that depend on maxnc 378 maxnc = maxnc + 50 ! this is not optimal 379 ! but the arrays could be shrunk 380 ! in routine ordern 381 call re_alloc(listc,1,maxnc,1,no_cl,name='listc') 382 call re_alloc(c,1,maxnc,1,no_cl,1,nspin,name='c') 383 endif 384 385 listc(nm+1,mul) = index ! LWF index 386 nct = nct + 1 ! number of coefficients so far 387 388 ! Find out structure of F and Ft matrices 389 ! F = S*C 390 ! find orbitals (nu, global index) which interact with mu 391 392 if (ncT2P(mul).gt.0) then 393 ! mul is one of 1:no_l 394 mull = ncT2P(mul) 395 indmu = listhptr(mull) 396 do imu = 1,numh(mull) 397 nu = listh(indmu+imu) 398 if (ilist(nu) .eq. 0) then ! Avoid overcounting 399 ilist(nu) = 1 400 numft(nu) = numft(nu) + 1 401 nf = nf + 1 402 endif 403 enddo 404 endif 405 ! At the end of this process, numft(nu) will contain 406 ! the number of locally managed orbitals (indexes 1:no_l) 407 ! that interact with nu. nf will be the total number of such interactions. 408 ! It is somehow the "sparsity" of the transpose of the S and H matrices 409 ! But note that numft is deallocated here, and only nftmax = max(numft(:)) 410 ! is retained. 411 412 ! Assign random guess for orbitals in atom ia if iopt = 0 413 ! Note: only those on atom ia, to avoid duplication of work 414 if (iopt .eq. 0) then 415 if (ja .eq. ia) then 416 call initguess(ia,iorb,isa(ia),nelectr, 417 . cg,cgval,Node) 418 do is = 1,nspin 419 c(nm+1,mul,is) = cg 420 enddo 421 snor = snor + cg**2 422 endif 423 endif 424 425 endif ! if iorb is part of our c-list 426 427 enddo loop_orbs_on_neighbor ! iorb 428 enddo loop_neighbor_atoms ! ja 429 430 ! Normalize LWF's if iopt = 0 ............................................. 431 ! (normalize to one if functions are expected to be occupied, 432 ! 0.5 if half occupied 433 ! 0.1 if empty) 434 ! 435 if (iopt .eq. 0) then 436 fact(1:2) = 1.0d0 437 if (ioptlwf .eq. 1) then 438 if (indexb.eq.indexi) then 439 if (nspin.eq.1) then 440 if (2*(nelectr/2) .eq. nelectr) fact(1:2) = sqrt(0.1) 441 if (2*(nelectr/2) .ne. nelectr) fact(1:2) = sqrt(0.5) 442 else 443 if (indexb.gt.nelectr1) fact(1) = sqrt(0.1) 444 if (indexb.gt.nelectr2) fact(2) = sqrt(0.1) 445 endif 446 endif 447 endif 448 449 do is = 1,nspin 450 do mu = lasto(ia-1)+1, lasto(ia) 451 mul = ncG2L(mu) 452 if (mul.gt.0) then ! mu is in our c-list 453 do in = 1,numc(mul) 454 if (listc(in,mul) .eq. index) then 455 !AG: can put spin loop here 456 c(in,mul,is) = c(in,mul,is) * fact(is)/sqrt(snor) 457 endif 458 enddo 459 endif 460 enddo 461 enddo 462 endif 463 464 nctmax = max ( nctmax , nct ) 465 nfmax = max ( nfmax , nf ) 466 467 enddo loop_lwfs_on_ia 468 enddo loop_ia 469 470 do mul = 1,no_cl 471 ncmax = max ( ncmax , numc(mul) ) 472 enddo 473 474 do mu = 1,no_u 475 nftmax = max ( nftmax , numft(mu) ) 476 enddo 477 478 nbands = index 479 480 if (index .gt. no_u) then 481 if (Node.eq.0) then 482 write(6,*) 'cspa: Number of LWFs larger than basis set size' 483 write(6,*) ' Increase basis set, or use less LWFs' 484 endif 485 call die() 486 endif 487 488 if ((ioptlwf .eq. 2) .and. (nbands .ne. nqtot/2)) then 489 if (Node.eq.0) then 490 write(6,*) 'cspa: Number of LWFs incorrectly calculated' 491 write(6,*) ' Something went wrong in generating the' 492 write(6,*) ' LWFs for the Ordejon-Mauri functional' 493 endif 494 call die() 495 endif 496 497 498C Find out sparse dimensions of Hij 499C loop over the localized wave funcions (centered at atoms).................. 500 501C Maximum interacion range between LWF's 502 r = 2.0 * rcoor + rh 503 504 if (2*r .ge. rmax) then 505 nhijmax = nbands 506 else 507 call mneighb(cell,r,natoms,xa,0,0,nna) !! ,nnmax) 508 509 index = 0 510 do ia = 1,natoms 511 512C Look for neighbours of atom ia within maximum interaction range 513 call mneighb(cell,r,natoms,xa,ia,0,nna) !! ,nnmax) 514 515 nhij = 0 516C Loop over the neighbors of ia within rcoor 517 do jj = 1,nna 518 ja = jan(jj) 519 alist(ja) = 0 520 enddo 521 do jj = 1,nna 522 ja = jan(jj) 523 if (alist(ja) .eq. 1) goto 20 524 alist(ja) = 1 525 526C determine how many LWF's centered in ja, depending on the atomic species 527C (or the number of electrons) 528 nelectr = qa (ja) + tiny 529 if (ioptlwf .eq. 1) then 530 indexj = ( ( nelectr + 2 ) / 2 ) 531 nelectr2 = nelectr/2 532 nelectr1 = nelectr - nelectr2 533 else if (ioptlwf .eq. 2) then 534 !AG: Is this right? 535 if ( (nelectr/2)*2 .ne. nelectr) then 536 call 537 $ die("Check Ordejon-Mauri Func. indexj assignment in cspa") 538 endif 539 indexj = ( ( nelectr / 2 ) ) 540 else 541 call die('cspa: Wrong functional option in cspa') 542 endif 543 544 nhij = nhij + indexj 54520 continue 546 enddo 547 do jj = 1,nna 548 ja = jan(jj) 549 alist(ja) = 0 550 enddo 551 nhijmax = max ( nhijmax , nhij ) 552 enddo ! ia 553 endif 554 555C Deallocate local memory 556 call reset_neighbour_arrays( ) 557 call de_alloc( lneeded, name='lneeded' ) 558 call de_alloc( numft, name='numft' ) 559 call de_alloc( ilist, name='ilist' ) 560 call de_alloc( alist, name='alist' ) 561 call de_alloc( indexloc, name='indexloc' ) 562 563 CONTAINS 564 565 !--------------------------------------------- 566 subroutine get_number_of_lwfs_on_atom() 567 568 ! Determine how many LWF's depending on the atomic species 569 ! (or the number of electrons) 570 nelectr = qa(ia) + tiny 571 if (abs(nelectr - qa(ia) + tiny) .gt. 1e-3) then 572 if (Node.eq.0) then 573 write(6,*) 'cspa: Wrong atomic charge for atom ',ia 574 write(6,*) ' qa = ',qa(ia),' must be an integer' 575 endif 576 call die() 577 endif 578 if (ioptlwf .eq. 1) then 579 indexi = ( ( nelectr + 2 ) / 2 ) 580 nelectr2 = nelectr/2 581 nelectr1 = nelectr - nelectr2 582 else if (ioptlwf .eq. 2) then 583 if ( (nelectr/2)*2 .ne. nelectr) then 584c EA: Instead of dying if there is any atom of odd species, use 585c an alternating scheme for assigning 1/2 more or 1/2 less, in 586c strict order of appearance. The user controls where the odd LWFs 587c go by defining the order of atoms in the AtomicSpeciesAndAtomicCoor... block 588c OLD------ if (Node.eq.0) then 589c write(6,*) 'cspa: Wrong Order-N functional option in ', 590c . 'cspa.' 591c write(6,*) ' You can only use the functional of' 592c write(6,*) ' Ordejon-Mauri for atoms with an even' 593c write(6,*) ' number of electrons.' 594c endif 595c OLD------ call die() 596c give one-extra/one-less LWF to odd species in turn with flag secondodd 597 if (secondodd) then 598 indexi = ( ( nelectr - 1 ) / 2 ) 599 secondodd = .false. 600 else 601 indexi = ( ( nelectr + 1 ) / 2 ) 602 secondodd = .true. 603 endif 604 else 605 indexi = ( ( nelectr ) / 2 ) 606 endif 607 else 608 call die('cspa: Wrong functional option in cspa') 609 endif 610 end subroutine get_number_of_lwfs_on_atom 611 612 end subroutine cspa 613 614 subroutine initguess(ia,iorb,is,ne,cg,cgval,Node) 615C ***************************************************************************** 616C Routine to assign an initial guess for an atomic orbital iorb in a localized 617C wave function centered on atom ia. 618C Assigns a random guess if the orbital belongs to the first 'zeta' of the 619C atom (lowest energy shell of its angular momentum), and if the angular 620C momentum is populated in the free atom. Otherwise, sets coefficient to cero. 621C 622C Written by P.Ordejon, November'96 623C lmax, lmaxs and nzls erased from the input by DSP, Aug. 1998. 624C ******************************* INPUT *************************************** 625C integer ia : Atom to which orbital belongs 626C integer mu : Orbital index within atom ia 627C integer is : Species index of atom ia 628C integer ne : Number of electrons of atom ia 629C real*8 cgval : Random value to set to cg if needed 630C integer Node : Node number 631C ***************************** OUTPUT **************************************** 632C real*8 cg : Initial guess for WF coefficient 633C ***************************************************************************** 634C The following functions must exist: 635C 636C INTEGER FUNCTION LOMAXFIS(IS) 637C Returns the maximum angular momentum of orbitals 638C Input: 639C INTEGER IS : Species index 640C 641C INTEGER FUNCTION NZTFL(IS,L) 642C Returns the number of different basis functions with the 643C same angular momentum L. 644C Input: 645C INTEGER IS : Species index 646C INTEGER L : Angular momentum 647C 648C ***************************************************************************** 649 650 use precision 651 652 use atmfuncs, only : lomaxfis, nztfl 653 use sys, only : die 654 655 implicit none 656 657 integer 658 . ia,iorb,is,ne,Node 659 660 real(dp) :: 661 . cg, cgval 662 663C Internal variables ......................................................... 664 integer 665 . index,iz,l,lmaxp,m 666 667C Initialize cg to cero 668 cg = 0.0d0 669 670C Find out angular momentum of orbital iorb 671 672 index = 0 673 do l = 0,lomaxfis(is) 674 do iz = 1,nztfl(is,l) 675 do m = -l,l 676 index = index + 1 677 enddo 678 if (index .ge. iorb) goto 10 679 enddo 680 enddo 681 call die('cspa: Error in orbital indexing in initguess') 68210 continue 683 684C Return if orbital is not the first zeta of its shell 685 if (iz .ne. 1) return 686 687C Assign initial guess. 688C If 2 or less electrons, populate lowest s orbital 689C If 8 or less electrons, populate lowest s and p orbitals 690C If 18 or less electrons, populate lowest s, p and d orbitals 691C If 32 or less electrons, populate lowest s, p, d and f orbitals 692 693 lmaxp = 0 694 if (ne .le. 32) lmaxp = 3 695 if (ne .le. 18) lmaxp = 2 696 if (ne .le. 8) lmaxp = 1 697 if (ne .le. 2) lmaxp = 0 698 if (ne .gt. 32) then 699 if (Node.eq.0) then 700 write(6,*) 'cspa: Cannot build initial guess in initguess.' 701 write(6,*) ' Reason: Too many electrons for this routine' 702 endif 703 call die() 704 endif 705 706 if (lmaxp .gt. lomaxfis(is)) then 707 if (Node.eq.0) then 708 write(6,*) 'cspa: Cannot build initial guess in initguess.' 709 write(6,*) ' Reason: Max. angular moment for atom ',ia, 710 . ' is not large enough' 711 endif 712 call die() 713 endif 714 715 if (ne .gt. 32) then 716 if (Node.eq.0) then 717 write(6,*) 'cspa: Cannot build initial guess in initguess.' 718 write(6,*) ' Too many valence electrons in atom ',ia 719 endif 720 call die() 721 endif 722 723 if (l .le. lmaxp) then 724 cg = cgval 725 endif 726 727 end subroutine initguess 728