CONGEN 5.0
Configuration generation for SCATCI
Loading...
Searching...
No Matches
congen_distribution.f90
Go to the documentation of this file.
1! Copyright 2019
2!
3! For a comprehensive list of the developers that contributed to these codes
4! see the UK-AMOR website.
5!
6! This file is part of UKRmol-in (UKRmol+ suite).
7!
8! UKRmol-in is free software: you can redistribute it and/or modify
9! it under the terms of the GNU General Public License as published by
10! the Free Software Foundation, either version 3 of the License, or
11! (at your option) any later version.
12!
13! UKRmol-in is distributed in the hope that it will be useful,
14! but WITHOUT ANY WARRANTY; without even the implied warranty of
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16! GNU General Public License for more details.
17!
18! You should have received a copy of the GNU General Public License
19! along with UKRmol-in (in source/COPYING). Alternatively, you can also visit
20! <https://www.gnu.org/licenses/>.
21!> \brief Distribute electrons to available orbitals
22!>
23module congen_distribution
24
25 implicit none
26
27 ! module entry point
28 public distrb
29
30 ! internal support routines
31 private assign
32 private cplea, cplem
33 private cgcoef
34 private getsa, getsm, getso
35 private packdet
36 private print1, print2, print3, print4, print5
37 private state
38 private wfn
39
40contains
41
42 !> \brief Assign quantum numbers to real shells.
43 !>
44 !> \param nshl Number of shells.
45 !> \param ndist Number of distributions generated from a given set of shells.
46 !> \param nused Number of slots used for finally assigned shells.
47 !> \param refcon ?
48 !> \param excon ?
49 !> \param qnstor Finally assigned quantum numbers.
50 !> \param nx Dimension of x overlayed with \c qnstor.
51 !> \param nftw Text output file unit.
52 !>
53 subroutine assign (nshl, ndist, nused, refcon, excon, qnstor, nx, nftw)
54
55 ! Reminder of some used global variables:
56 ! pqn(1,nshl) 0 for pseudoshell, sequence number for real shell
57 ! pqn(2,nshl) starting index for pseudoshell
58 ! pqn(3,nshl) ending index for pseudoshell
59 ! occshl(nshl) occupation
60 ! sshl(nshl) sym-values per shell
61 ! nst(m,n) pointer to shell index and order count
62 ! nshsym(nsym) counter for number of shells in s symmetry
63 ! qnshlr(nshl) work area stores index of shells
64
65 use congen_data, only : ntcon, test, nrcon, nobt, nst, nshsym, nobi, ift, occshl, pqn => pqnst, sshl => sshlst, qnshlr
66
67 integer :: ndist, nftw, nshl, nused, nx
68 integer, dimension(nobt) :: excon
69 integer, dimension(*) :: qnstor
70 integer, dimension(nobt,max(ntcon,1)) :: refcon ! JMC changing the 2nd dimension from 2
71 intent (in) nftw, nshl, nx, refcon
72 intent (out) nused, qnstor
73 intent (inout) excon, ndist
74
75 integer, save :: allow, i, ic, it, ita, itest, j, ks, kss, mrun, mt, nav, navm1, ndrop, nrep, nshrun, nt
76
77 ! workspace limits
78 if (ift /= 0) then
79 nav = nx
80 navm1 = nav - 1
81 ift = 0
82 end if
83
84 ! categorize shells by symmetry (into NST array)
85
86 mrun = 0 ! 1-based index of current symmetry being processed (-> number of symmetries processed)
87 nshrun = 0 ! 1-based index of current shell being processed (-> number of shells proccessed)
88
89 do while (nshrun /= nshl) ! process all NSHL shells
90 mrun = mrun + 1
91 ita = 0 ! index of shell within symmetry (-> number of shells within symmetry)
92 do i = 1, nshl
93 if (mrun == sshl(i)) then
94 ita = ita + 1
95 nst(mrun,ita) = i
96 end if
97 end do
98 nshrun = nshrun + ita
99 end do
100
101 ! initializations
102
103 nshsym(1:mrun) = 0 ! number of shells processed so far per symmetry
104 ndist = 0 ! number of distributions generated (since the start of the program)
105 ndrop = 0 ! number of distributed shells stored (since the start of the program)
106 ks = 0 ! 1-based index of the current shell
107
108 ! ???
109
110 100 ks = ks + 1 ! next shell, please!
111 mt = sshl(ks) ! get symmetry number of the current shell
112 it = nshsym(mt) ! get number of shells with the same symmetry processed before the current one
113 nshsym(mt) = it + 1 ! update number of processed shells with that symmetry
114 qnshlr(ks) = pqn(2,ks) ! starting orbital of the set of orbitals this shell belongs to
115
116 ! skip to 200 if any of the previous orbitals for this symmetry violate rules below
117 120 do i = 1, it
118 kss = nst(mt,i)
119 if (qnshlr(kss) == qnshlr(ks)) go to 200
120 if (occshl(kss) == occshl(ks) .and. &
121 qnshlr(kss) >= qnshlr(ks) .and. &
122 qnshlr(ks) >= pqn(2,kss) .and. &
123 qnshlr(kss) <= pqn(3,ks)) go to 200
124 end do
125
126 if (ks < nshl) go to 100
127
128 ! take into account possible constraints
129 if (ntcon /= 0) then
130 excon(1:nobt) = 0
131 do i = 1, nshl
132 j = nobi(sshl(i)) + qnshlr(i)
133 excon(j) = occshl(i)
134 end do
135
136 allow = 0
137 itest = 0
138
139 do ic = 1, ntcon
140 nrep = 0
141 do i = 1, nobt
142 nt = excon(i) - refcon(i,ic)
143 nrep = abs(nt) + nrep
144 end do
145 if (test(ic) /= 1) then
146 if (nrep <= nrcon(ic) + nrcon(ic)) go to 200
147 else
148 if (nrep <= nrcon(ic) + nrcon(ic)) allow = 1
149 end if
150 itest = itest + test(ic)
151 end do
152
153 if (allow == 0 .and. itest /= 0) go to 200
154 end if
155
156 ! allowed assignment is to be stored - check that there is enough space
157 if (ndrop + nshl >= navm1) then
158 write(nftw,'("1",31("*"),/," ",31("*"),"STORAGE OVERFLOW IN ASSIGN:",I8, " WORDS AVAILABLE")') nx
159 stop 70
160 end if
161
162 ! store the assignment
163 ndist = ndist + 1
164 do i = 1, nshl
165 ndrop = ndrop + 1
166 qnstor(ndrop) = qnshlr(i)
167 end do
168
169 ! ascend in assignment loops
170 200 do while (ks > 0)
171 mt = sshl(ks) ! get symmetry number of the current shell
172 it = nshsym(mt) - 1 ! get number of shells with the same symmetry processed before the current one
173 qnshlr(ks) = qnshlr(ks) + 1
174 if (qnshlr(ks) <= pqn(3,ks)) go to 120
175 nshsym(mt) = it
176 ks = ks - 1
177 end do
178
179 nused = ndrop
180
181 end subroutine assign
182
183
184 !> \brief Clebsch-Gordan coefficients.
185 !>
186 !> Evaluates all Clebsch-Gordan coefficients for given \f$ j_1, j_2, j_2 \f$ and \f$ m_3 \f$.
187 !>
188 !> \param j1 Angular momentum 1.
189 !> \param j2 Angular momentum 2.
190 !> \param j3 Angular momentum 3.
191 !> \param m3 Angular momentum projection 3.
192 !> \param n On return, number of (\c m1, \c m2) pairs giving non-zero Clebsch-Gordan coefficients.
193 !> \param ms On return, pairs of projections \c m1, \c m2 for angular momenta j1 and j2.
194 !> \param c On return, array of resulting \c n Clebsch-Gordan coefficients.
195 !> \param intpfg Requested text output intensity.
196 !>
197 subroutine cgcoef (j1, j2, j3, m3, n, ms, c, intpfg)
198
199 use precisn, only : wp
200 use consts, only : xzero, xone
201 use congen_data, only : binom, ind, nftw, thresh1
202
203 integer :: intpfg, j1, j2, j3, m3, n
204 real(kind=wp), dimension(*) :: c
205 integer, dimension(2,*) :: ms
206 intent (in) intpfg, j1, j2, j3, m3
207 intent (inout) c, ms, n
208
209 real(kind=wp) :: a, b, t
210 integer :: i, i1, i2, ii, jj, js, lb, lb1, lb2, lbh, lbl, m, m1, m2
211 integer, dimension(3) :: j, k, l
212
213 js = (j1 + j2 + j3 - 1) / 2
214 j(1) = j1
215 j(2) = j2
216 j(3) = j3
217 k(1) = js - j2
218 k(2) = js - j3
219 k(3) = js - j1
220 n = 0
221 m = m3
222 if (j3 - 1 < abs(m - 1)) return
223 if (any(k < 0)) return
224 a = xone / (binom(ind(js+1)+k(2)) * binom(ind(j3)+k(1)))
225 l(3) = (j3 + m3 - 2) / 2
226 m1 = j1
227
228 do jj = 1, j1
229 m2 = m3 - m1 + 1
230 if (abs(m2 - 1) <= j2 - 1) then
231 n = n + 1
232 ms(1,n) = m1
233 ms(2,n) = m2
234 l(1) = (j1 - m1) / 2
235 l(2) = (j2 + m2 - 2) / 2
236 b = a
237 do ii = 1, 3
238 b = binom(ind(j(ii)) + k(ii)) / binom(ind(j(ii)) + l(ii)) * b
239 end do
240 b = sqrt(b)
241 i1 = max(l(1) - k(1), l(2) - k(3), 0)
242 i2 = min(l(1), l(2), k(2))
243 t = xzero
244 if (i2 >= i1) then
245 lbl = ind(k(2) + 1) + i1
246 lbh = ind(k(2) + 1) + i2
247 lb1 = ind(k(1) + 1) + l(1) - i1
248 lb2 = ind(k(3) + 1) + l(2) - i1
249 do lb = lbl, lbh
250 t = binom(lb) * binom(lb1) * binom(lb2) - t
251 lb1 = lb1 - 1
252 lb2 = lb2 - 1
253 end do
254 end if
255 c(n) = b * t * (-xone)**i2
256 if (abs(c(n)) <= thresh1) n = n - 1
257 end if
258 m1 = m1 - 2
259 end do
260
261 if (intpfg /= 0) then
262 write(nftw,'(" CGCOEF : CLEBSCH-GORDAN COEFFICIENTS FOR")')
263 write(nftw,'(" J1 =",I4," J2 =",I4," J3 =",I4," M3 =",I4,/)') j1, j2, j3, m3
264 write(nftw, '(/(E25.15,2I5))') (c(i), ms(1,i), ms(2,i), i=1,n)
265 end if
266
267 end subroutine cgcoef
268
269
270 !> \brief Form \c ne electrons in a shell coupled to (l,is).
271 !>
272 !> This version is for linear molecules (Alchemy).
273 !>
274 !> \param ne Number of electrons.
275 !> \param l Lambda of shell.
276 !> \param is Spin of coupled shell (is=2*s+1).
277 !> \param isz 2*sz+1
278 !> \param m Projection of lambda of coupled shell.
279 !> \param nc Number of determinants required for shell.
280 !> \param c Coefficient of determinant,
281 !> \param iso Spin-orbitals for determinants according to the following table
282 !> \verbatim
283 !> 0 sa 0 l+a
284 !> 1 sb 1 l+b
285 !> 2 l-a
286 !> 3 l-b
287 !> \endverbatim
288 !>
289 subroutine getsa (ne, l, is, isz, m, nc, c, iso)
290
291 use precisn, only : wp
292 use consts, only : one => xone
293 use congen_data, only : root2
294
295 integer :: is, isz, l, m, nc, ne
296 real(kind=wp), dimension(*) :: c
297 integer, dimension(*) :: iso
298 intent (in) is, isz, l, m, ne
299 intent (out) c, nc
300 intent (inout) iso
301
302 nc = 1
303 c(1) = one
304 if (ne == 0) return
305 if (ne /= 4) then
306 if (ne < 2) then
307 iso(1) = 1 - isz / 2
308 if (l == 0) return
309 if (m < 0) iso(1) = iso(1) + 2
310 return
311 else if (ne == 2) then
312 if (l == 0) then
313 iso(1) = 0
314 iso(2) = 1
315 return
316 end if
317 if (isz + m /= 1) then
318 if (m /= 0) then
319 iso(1) = 1 - (l + l) / m
320 iso(2) = iso(1) + 1
321 return
322 end if
323 iso(1) = (3 - isz) / 4
324 iso(2) = iso(1) + 2
325 return
326 end if
327 nc = 2
328 iso(1) = 0
329 iso(2) = 3
330 iso(3) = 1
331 iso(4) = 2
332 c(1) = root2
333 c(2) = root2
334 if (is == 1) c(2) = -root2
335 return
336 else
337 if (m < 0) then
338 iso(1) = 1 - isz / 2
339 iso(2) = 2
340 iso(3) = 3
341 return
342 end if
343 iso(1) = 0
344 iso(2) = 1
345 iso(3) = 3 - isz / 2
346 return
347 end if
348 end if
349
350 iso(3) = 2
351 iso(4) = 3
352 iso(1) = 0
353 iso(2) = 1
354
355 end subroutine getsa
356
357
358 !> \brief Form \c ne electrons in a shell coupled to (l,is).
359 !>
360 !> This version is for non-linear molecules (Molecule).
361 !>
362 !> \param ne Number of electrons
363 !> \param isz 2*sz+1
364 !> \param nc Number of determinants required for shell.
365 !> \param c Coefficient of determinant.
366 !> \param iso Spin-orbitals for determinants according to the following table:
367 !> \verbatim
368 !> 0 sa 0 l+a
369 !> 1 sb 1 l+b
370 !> 2 l-a
371 !> 3 l-b
372 !> \endverbatim
373 !>
374 subroutine getsm (ne, isz, nc, c, iso)
375
376 use precisn, only : wp
377 use consts, only : one => xone
378
379 integer :: isz, nc, ne
380 real(kind=wp), dimension(*) :: c
381 integer, dimension(*) :: iso
382 intent (in) isz, ne
383 intent (out) c, iso, nc
384
385 nc = 1
386 c(1) = one
387 iso(1) = 0
388 iso(2) = 1
389 if (ne == 1) iso(1) = 1 - isz / 2
390
391 ! So the results will be:
392 !
393 ! isz ne iso(1) iso(2)
394 ! ---------------------------------
395 ! any /= 1 0 1
396 ! 1 == 1 1 1
397 ! 2 == 1 1 1
398 ! 3 == 1 0 1
399 ! 4 == 1 -1 1
400
401 end subroutine getsm
402
403
404 !> \brief ?
405 !>
406 !> \param ns ?
407 !> \param intpfg Print level flag.
408 !> \param nti ?
409 !> \param iqns ?
410 !> \param ci ?
411 !> \param nd ?
412 !> \param id ?
413 !> \param cd ?
414 !> \param last ?
415 !>
416 subroutine getso (ns, intpfg, nti, iqns, ci, nd, id, cd, last)
417
418 use precisn, only : wp
419 use congen_data, only : nes => occshl, ms => mshl, iqn => qnshl, ne => nnlecg, symtyp, nftw
420
421 integer :: intpfg, last, nd, ns, nti
422 real(kind=wp), dimension(*) :: cd, ci
423 integer, dimension(*) :: id
424 integer, dimension(2,ns,*) :: iqns
425 intent (in) ci, intpfg, iqns, last, ns, nti
426 intent (inout) cd, id, nd
427
428 real(kind=wp), dimension(100) :: c ! JMC change dimension to ns*(max nc=2) ??? (an overestimate, for safety).
429 real(kind=wp), dimension(ns) :: cs ! JMC changing the dimension from 50
430 integer :: i, ie1, ie2, is, isz, iti, kc, ke, kso, ld, ld1, ld2, ml, nc
431 integer, dimension(200) :: iso ! JMC change dimension to (max nc=2)*sum(nes(i),i=1,ns) ??? (an overestimate, for safety).
432 integer, dimension(150) :: jso ! JMC change dimension to max(ne, sum(nes(i),i=1,ns)) ???
433 integer, dimension(ns+1) :: lc, lso ! JMC changing the dimension from 51. Think LSO could be (NS) not (NS+1)...
434 integer, dimension(ns) :: lcs, lsos ! JMC changing the dimension from 50
435 real(kind=wp) :: t
436
437 if (intpfg /= 0) then
438 write(nftw,'(" GETSO : NTI =",I10,/," IQNS :",/)') nti
439 do iti = 1, nti
440 write(nftw,'(20I5)') (iqns(1,is,iti), is=1,ns)
441 write(nftw,'(20I5)') (iqns(2,is,iti), is=1,ns)
442 end do
443 end if
444
445 nd = 0
446 ld = 1
447
448 do iti = 1, nti
449
450 kc = 1
451 kso = 1
452
453 do is = 1, ns
454 isz = iqns(1,is,iti)
455 if (symtyp <= 1) then
456 ml = iqns(2,is,iti)
457 call getsa (nes(is), ms(is), iqn(1,is), isz, ml, nc, c(kc), iso(kso))
458 else
459 call getsm (nes(is), isz, nc, c(kc), iso(kso))
460 end if
461 lc(is) = kc
462 lso(is) = kso
463 kc = kc + nc
464 kso = kso + nc * nes(is)
465 end do
466
467 lc(ns + 1) = kc
468 is = 1
469 t = ci(iti)
470 ie2 = 0
471
472 300 cs(is) = t
473 lsos(is) = lso(is)
474 lcs(is) = lc(is)
475
476 if (nes(is) == 0) go to 415
477
478 400 ie1 = ie2 + 1
479 ie2 = ie2 + nes(is)
480 kso = lsos(is)
481 do ke = ie1, ie2
482 jso(ke) = iso(kso)
483 kso = kso + 1
484 end do
485
486 415 lsos(is) = kso
487 t = cs(is) * c(lcs(is))
488 lcs(is) = lcs(is) + 1
489 is = is + 1
490 if (is <= ns) go to 300
491
492 nd = nd + 1
493 if (nd > last) then
494 write(nftw,'("0storage overflow")')
495 nd = 0
496 return
497 end if
498
499 do ke = 1, ne
500 id(ld) = jso(ke)
501 ld = ld + 1
502 end do
503
504 cd(nd) = t
505
506 do while (is > 1)
507 is = is - 1
508 ie2 = ie2 - nes(is)
509 if (lcs(is) < lc(is + 1)) go to 400
510 end do
511
512 end do
513
514 if (intpfg /= 0) then
515 write(nftw,'(" GETSO : ND =",I6,/," CD, ID :",/)') nd
516 ld2 = 0
517 do i = 1, nd
518 ld1 = ld2 + 1
519 ld2 = ld2 + ne
520 write(nftw,'(E25.15,20I5)') cd(i), (id(ld), ld=ld1,ld2)
521 end do
522 end if
523
524 end subroutine getso
525
526
527 !> \brief Loop through (and fill) all allowed couplings for a given electron distribution into shells.
528 !>
529 !> \verbatim
530 !> CUPSET IS THE ENTRY TO SET LOCATIONS OF ARRAYS
531 !> MSHL(NSHL) SYMMETRY NUMBER FROM ZERO TO N-1 (MVALUE)
532 !> QNSHL(3,2*NSHL-1)
533 !> 1 -- MULT / 2 -- SYMMETRY / 3 -- +- (NOT USED)
534 !> CUP(3,NSHL-1)
535 !> QNTOT(3) TOTAL QN'S
536 !> SPNMIN(NSHL-1) TEMP STORAGE FOR LOWEST SPIN COUPLING
537 !> X(NX) WORK AREA (R*8)
538 !> NSHL NUMBER OF TRUE SHELLS OCCUPIED
539 !> NSTATE NUMBER OF COUPLINGS (COMPUTED)
540 !> NTYPE PROTO-TYPE NUMBER (INPUT)
541 !> NDIST NUMBER OF PQN ASSIGNMENTS
542 !> NCSF RUNNING CSF NUMBER
543 !> SYMTYP GE 2 FOR MOLECULE
544 !> CONPF PRINT FLAG
545 !>
546 !> GUSHELL(NSHL) GU VALUE FOR EACH SHELL
547 !> SHLMX(NSHL) MAX OCCUPATION FOR A SHELL
548 !> OCCSHL(NSHL) COMPRESSED SHELL OCC'S / ALL ZEROS DELETED
549 !> AND PSEUDO SHELLS EXPANDED
550 !> QNSHL(3,NSHL) FIRST INDEX = 1 MULTIPLICITY FOR EACH COUPLING
551 !> FIRST INDEX = 2 M-VALUE FOR EACH COUPLING
552 !> FIRST INDEX = 3 +1 FOR (S+) -1 FOR (S-)
553 !> OTHERWIZE ZERO
554 !> CUP(3,NSHL) COUPLING SCEME IN ORDER STATE A TO STATE B
555 !> GIVES STATE C
556 !> QNTOT(3) INPUT STATE TO BE SEARCHED FOR. ORDER: MULTI-
557 !> PLICITY,ANG.MOM., PLUS(+1) OR MINUS(-1)
558 !> GUTOT G(+1) OR U(-1) FOR INPUT STATE
559 !> KSLIM(2,NSHL) FIRST INDEX = 1 IS NUMBER OF STATES
560 !> ALREDY LOOPED OVER IN A SHELL
561 !> FIRST INDEX = 2 IS THE NUMBER OF STATES FOR
562 !> A SHELL
563 !> MCLIM(2,NSHL) FIRST INDEX = 1 IS THE NUMBER OF STATES
564 !> ALREADY LOOPED OVER IN A COUPLING
565 !> FIRST INDEX = 2 IS THE MAXIMUM NUMBER OF STATES
566 !> IN A COUPLING
567 !> QNTMP(2,3,NSHL) FIRST INDEX POINTS ON + OR
568 !> QNTMP(2,3,NSHL) FIRST INDEX POINTS ON +- AND M(2 AND 1 RESP.)
569 !> SECOND INDEX POINTS ON POSSIBLE COUPLINGS
570 !> (M+M,M-M,2*M,S+,S-) SAVE AREA DOWN THE LOOPS
571 !> SPNMIN(NSHL) MIN MULTIPLICITY FOR A COUPLING
572 !> \endverbatim
573 !>
574 subroutine cplea (nncsf, nadel, x, nftw)
575
576 use iso_c_binding, only : c_loc, c_f_pointer
577 use precisn, only : wp
578 use congen_data, only : lg, next, nx, icdi, iexcon, indi, inodi, iqnstr, irfcon, ndel, confpf, ndist, nshl, &
579 nstate, occshl, mshl, gushl, cup, qnshl, spnmin, kslim, mclim, shlmx1, qntot, &
580 qntar, gutot, nndel, cdimx, ndimx, nodimx
581
582 integer :: nadel, nftw, nncsf
583 real(wp), dimension(*), target :: x
584 integer, dimension(:), pointer :: int_ptr, irfcon_ptr, iexcon_ptr, iqnstr_ptr, nodi_ptr, ndi_ptr, ndel_ptr, nnext_ptr
585 integer, dimension(2,3,lg) :: qntmp
586
587 integer :: gutry, i, idumm, iidis1, iidist, iocc, kc, kclim, kcs, kcsp, ks, &
588 m, mtry, nav, nnext, nused, shl1, shl2, spntry, z1=1
589
590 real(wp), pointer :: x_ptr
591
592 nstate = 0 ! number of states which may be formed with current distribution into shells
593 iidis1 = 0 ! number of the current shell distribution
594
595 ! perform gross checks on symmetries allowed for dist
596
597 if (gutot /= 0) then
598 gutry = 1
599 do i = 1, nshl
600 if (gushl(i) > 0) cycle
601 if (iand(occshl(i), z1) /= 0) gutry = -gutry
602 end do
603 if (gutry /= gutot) return
604 end if
605 spntry = 1
606 mtry = 0
607 do i = 1, nshl
608 m = mshl(i) + 1
609 if (gutot /= 0 ) m = m + m
610 iocc = occshl(i)
611 if (iocc > shlmx1(m) / 2) iocc = shlmx1(m) - iocc
612 spntry = spntry + iocc
613 mtry = mtry + iocc * mshl(i)
614 end do
615 if (spntry > qntot(1)) return
616 if (mtry - qntot(2) > 0 .or. iand(mtry - qntot(2), z1) /= 0) return
617
618 ! initialize third shell coupling
619
620 kclim = nshl - 1
621 do i = 1, kclim
622 qntmp(1,3,i) = 0
623 qntmp(2,3,i) = -1
624 end do
625
626 ! initialize shell quantum numbers and loop limits for shell coupling
627
628 spntry = 1
629 mtry = 0
630 do i = 1, nshl
631 m = mshl(i) + 1
632 if (gutot /= 0) m = m + m
633 qnshl(1,i) = 1
634 qnshl(2,i) = 0
635 qnshl(3,i) = 1
636 kslim(2,i) = 1
637 if (occshl(i) == 0 .or. occshl(i) == shlmx1(m)) cycle
638 spntry = spntry + 1
639 mtry = mtry + mshl(i)
640 qnshl(1,i) = 2
641 qnshl(2,i) = mshl(i)
642 qnshl(3,i) = 0
643 if (mshl(i) == 0) qnshl(3,i) = 1
644 if (occshl(i) /= 2) cycle
645 kslim(2,i) = 3
646 spntry = spntry - 1
647 mtry = mtry - mshl(i)
648 end do
649
650 ! begin to decend into loops
651
652 ks = 1
653 100 kslim(1,ks) = 1
654 if (kslim(2,ks) == 1) go to 120
655 qnshl(1,ks) = 3
656 qnshl(2,ks) = 0
657 qnshl(3,ks) = -1
658 spntry = spntry + 2
659 120 ks = ks + 1
660 if (ks < nshl) go to 100
661
662 if (spntry < qntot(1) .or. mtry < qntot(2)) go to 600
663 kc = 1
664 kcs = nshl + 1
665 if (kclim /= 0) go to 300
666 if (qnshl(1,1) /= qntot(1) .or. qnshl(2,1) /= qntot(2)) go to 500
667 if (qnshl(2,1) == 0 .and. qnshl(3,1) /= qntot(3)) go to 500
668 go to 420
669 300 shl1 = cup(1,kc)
670 shl2 = cup(2,kc)
671 qnshl(1,kcs) = qnshl(1,shl1) + qnshl(1,shl2) - 1
672 spnmin(kc) = abs(qnshl(1,shl1) - qnshl(1,shl2)) + 1
673 qntmp(1,1,kc) = qnshl(2,shl1) + qnshl(2,shl2)
674 qnshl(2,kcs) = qnshl(2,shl1) + qnshl(2,shl2)
675 qntmp(1,2,kc) = abs(qnshl(2,shl1) - qnshl(2,shl2))
676 qntmp(2,1,kc) = qnshl(3,shl1) * qnshl(3,shl2)
677 qnshl(3,kcs) = qnshl(3,shl1) * qnshl(3,shl2)
678 qntmp(2,2,kc) = 0
679 mclim(1,kc) = 1
680 mclim(2,kc) = 1
681 if (qntmp(1,1,kc) == qntmp(1,2,kc)) go to 320
682 mclim(2,kc) = 2
683 if (qntmp(1,2,kc) /= 0) go to 320
684 mclim(2,kc) = 3
685 qntmp(2,2,kc) = 1
686 320 kc = kc + 1
687 kcs = kcs + 1
688 if (kc < kclim) go to 300
689
690 ! test if the final coupling is in the range permitted
691
692 kc = kclim
693 kcs = kcs - 1
694 if (qnshl(1,kcs) < qntot(1) .or. spnmin(kc) > qntot(1)) go to 500
695 qnshl(1,kcs) = qntot(1)
696 if (qnshl(2,kcs) /= qntot(2) .and. qntmp(1,2,kc) /= qntot(2)) go to 500
697 qnshl(2,kcs) = qntot(2)
698 qnshl(3,kcs) = qntot(3)
699 if (qntmp(1,1,kc) /= 0) go to 420
700 qnshl(3,kcs) = qntmp(2,1,kc)
701 if (qnshl(3,kcs) /= qntot(3)) go to 500
702
703 ! for R-matrix calculations: reject couplings which do not preserve the target quantum numbers
704
705 420 if (qntar(1) >= 0) then
706 kcsp = kcs - 1
707 if (kcs == 3) kcsp = 1
708 if (qnshl(1,kcsp) /= qntar(1)) go to 500
709 if (qnshl(2,kcsp) /= qntar(2)) go to 500
710 if (qnshl(3,kcsp) /= qntar(3)) go to 500
711 end if
712 nstate = nstate + 1
713 if (nstate /= 1) go to 400
714
715 ! convert some real pointers to integer pointers for use in ASSIGN
716 x_ptr => x(irfcon) ; call c_f_pointer (c_loc(x_ptr), irfcon_ptr, (/1/))
717 x_ptr => x(iexcon) ; call c_f_pointer (c_loc(x_ptr), iexcon_ptr, (/1/))
718 x_ptr => x(iqnstr) ; call c_f_pointer (c_loc(x_ptr), iqnstr_ptr, (/1/))
719
720 ! Assign pqn value and print type and distrib data for allowed state
721 call assign (nshl, ndist, nused, irfcon_ptr, iexcon_ptr, iqnstr_ptr, nx, nftw)
722
723 if (ndist == 0) then
724 nstate = 0
725 return
726 end if
727
728 nav = nx - nused
729 nnext = next + nused
730 if (nndel /= 0) go to 405
731 call print1 (iqnstr_ptr)
732 400 if (confpf >= 10 .and. nndel == 0) call print2 (0)
733 405 continue
734
735 ! convert some real pointers to integer pointers for use in WFN
736 x_ptr => x(inodi) ; call c_f_pointer (c_loc(x_ptr), nodi_ptr, (/1/))
737 x_ptr => x(indi) ; call c_f_pointer (c_loc(x_ptr), ndi_ptr, (/1/))
738 x_ptr => x(ndel) ; call c_f_pointer (c_loc(x_ptr), ndel_ptr, (/1/))
739 x_ptr => x(iqnstr) ; call c_f_pointer (c_loc(x_ptr), iqnstr_ptr, (/1/))
740 x_ptr => x(nnext) ; call c_f_pointer (c_loc(x_ptr), nnext_ptr, (/1/))
741
742 call wfn (nncsf, nadel, iidist, iidis1, nodi_ptr, ndi_ptr, x(icdi), ndel_ptr, iqnstr_ptr, x(nnext), nnext_ptr, nav)
743
744 if (nndel > 0 .and. nadel > nndel .and. iidis1 == 0) return
745 if (nndel == 0) go to 500
746 if (iidist == 0) go to 500
747 if (nstate == 1) call print1 (iqnstr_ptr)
748 if (confpf >= 10 .and. iidis1 /= 0) call print2 (iidis1)
749
750 ! ascend in coupling tree - ascend in the shell to shell coupling loops
751
752 500 kc = kc - 1
753 kcs = kcs - 1
754 if (kc == 0) go to 600
755 mclim(1,kc) = mclim(1,kc) + 1
756 if (mclim(1,kc) <= mclim(2,kc)) go to 520
757 qnshl(1,kcs) = qnshl(1,kcs) - 2
758 if (qnshl(1,kcs) < spnmin(kc)) go to 500
759 mclim(1,kc) = 1
760 520 qnshl(2,kcs) = qntmp(1,mclim(1,kc),kc)
761 qnshl(3,kcs) = qntmp(2,mclim(1,kc),kc)
762 go to 320
763
764 ! ascend in the loops which couple shells to themselves
765
766 600 ks = ks - 1
767 if (ks == 0) go to 900
768 kslim(1,ks) = kslim(1,ks) + 1
769 if (kslim(1,ks) > kslim(2,ks)) go to 600
770 if (kslim(1,ks) == 3) go to 625
771
772 ! couple to 1(2*l) and adjust spntry and mtry
773
774 qnshl(1,ks) = 1
775 spntry = spntry - 2
776 qnshl(2,ks) = mshl(ks) + mshl(ks)
777 mtry = mtry + mshl(ks) + mshl(ks)
778 qnshl(3,ks) = 0
779 go to 120
780
781 ! couple to 1(s+) and adjust spntry and mtry
782
783 625 mtry = mtry - qnshl(2,ks)
784 qnshl(2,ks) = 0
785 qnshl(3,ks) = 1
786 go to 120
787
788 900 if (nndel /= 0 .and. iidis1 /= 0) call print3 (iidis1, 1)
789 if (nndel == 0 .and. nstate /= 0) call print3 (idumm, 2)
790
791 end subroutine cplea
792
793
794 !> \brief ?
795 !>
796 !> \param nncsf Number of CSFs (total).
797 !> \param nadel Number of processed configurations.
798 !> \param x Real workspace (or the free part of it).
799 !> \param nftw Text output unit.
800 !>
801 subroutine cplem (nncsf, nadel, x, nftw, noex)
802
803 use iso_c_binding, only : c_loc, c_f_pointer
804 use precisn, only : wp
805 use global_utils, only : mprod
806 use congen_data, only : lg, next, nx, icdi, iexcon, indi, inodi, iqnstr, irfcon, ndel, confpf, ndist, nshl, &
807 nstate, occshl, mshl, cup, qnshl, spnmin, shlmx1, qntot, qntar, nndel
808
809 ! Reminder of some used global variables:
810 !
811 ! occshl(nshl) compressed shell occ's / all zeros deleted and pseudo shells expanded
812 ! mshl(nshl) symmetry number from zero to n-1 (mvalue)
813 ! qnshl(3,2*nshl-1) 1 -- mult / 2 -- symmetry / 3 -- +- (not used)
814 ! cup(3,nshl-1)
815 ! qntot(3) total qn's
816 ! qntmp(2*nshl-1) temp storage for true molecule sym values -- zeros passed to bowen in qnshl
817 ! spnmin(nshl-1) temp storage for lowest spin coupling
818 ! x(nx) real work area
819 ! nshl number of true shells occupied
820 ! nstate number of couplings (computed)
821 ! ntype proto-type number (input)
822 ! ndist number of pqn assignments
823 ! ncsf running csf number
824 ! symtyp >= 2 for molecule
825 ! conpf print flag
826
827 integer :: nadel, nftw, nncsf, noex
828 real(wp), dimension(*), target :: x
829 integer, dimension(:), pointer :: int_ptr, nodi_ptr, ndi_ptr, ndel_ptr, irfcon_ptr, iexcon_ptr, iqnstr_ptr, nnext_ptr
830 real(wp), pointer :: x_ptr
831
832 integer :: i, idumm, iidis1, iidist, kc, kclim, kcs, kcsp, mtry, mu1, mu2, nav, nnext, nused, spntry
833 integer, dimension(2*lg) :: qntmp ! JMC increasing the dimension from LG to match 2nd dimension of QNSHL
834
835 ! initialize
836
837 nstate = 0 ! ?
838 mtry = 0 ! combination of IRRs of all shells
839 spntry = 1 ! maximal combination of spins of all shels
840 iidis1 = 0 ! ?
841
842 do i = 1, nshl ! loop over all shells in the current set of orbitals
843 qnshl(1,i) = 1 ! default to singlet (2S+1)
844 qnshl(2,i) = 0 ! default to maximal symmetry (IRR = 1)
845 qnshl(3,i) = 1 ! default to +1 inversion symmetry
846 qntmp(i) = 0
847 if (occshl(i) == shlmx1(mshl(i) + 1)) cycle ! fully occupied shells do not contribute anything valuable -> skip
848 qnshl(1,i) = 2 ! but shells not fully occupied are doublets
849 spntry = spntry + 1 ! compose maximal possible total spin
850 qntmp(i) = mshl(i) ! the symmetry number of the shell
851 mtry = mprod(mtry + 1, qntmp(i) + 1, 0, nftw) - 1 ! combine IRR of all shells so far
852 end do
853
854 if (spntry < qntot(1)) return ! jump out if impossible to compose the wanted total spin
855 if (mtry /= qntot(2)) return ! jump out if the combined IRR of all shells is not the wanted one
856 if (nshl == 1 .and. qnshl(1,1) /= qntot(1)) return ! jump out if the *only* shell violates spin requirement
857
858 ! complete M-values for coupling of shells to each other
859
860 kc = 0
861 kcs = nshl + 1
862
863 if (nshl == 1 .and. qnshl(1,1) == qntot(1)) go to 360 ! trivial combination -- one contributing shell
864
865 do i = 1, nshl - 1
866 qntmp(kcs) = mprod(qntmp(cup(1,i)) + 1, qntmp(cup(2,i)) + 1, 0, nftw) - 1
867 qnshl(2,kcs) = 0
868 qnshl(3,kcs) = 1
869 kcs = kcs + 1
870 end do
871
872 ! begin to descend into shell to shell coupling tree
873
874 kcs = nshl
875
876 300 do while (kc < nshl - 1)
877 kc = kc + 1
878 kcs = kcs + 1
879 mu1 = qnshl(1,cup(1,kc))
880 mu2 = qnshl(1,cup(2,kc))
881 qnshl(1,kcs) = mu1 + mu2 - 1
882 spnmin(kc) = abs(mu1 - mu2) + 1
883 end do
884
885 ! test for allowed state
886
887 if (qnshl(1,kcs) < qntot(1) .or. spnmin(kc) > qntot(1)) go to 500
888 qnshl(1,kcs) = qntot(1)
889
890 ! for R-matrix calculations: reject couplings which do not preserve the target quantum numbers
891
892 360 if (qntar(1) >= 0) then
893 kcsp = kcs - 1
894 if (kcs == 3) kcsp = 1
895 if (qnshl(1,kcsp) /= qntar(1)) go to 500
896 end if
897
898 nstate = nstate + 1
899
900 ! first iteration: assign pqn value and print type and distrib data for allowed state
901
902 if (nstate == 1) then
903 iexcon = 1
904
905 ! convert some real pointers to integer pointers for use in ASSIGN
906 x_ptr => x(irfcon) ; call c_f_pointer (c_loc(x_ptr), irfcon_ptr, (/1/))
907 x_ptr => x(iexcon) ; call c_f_pointer (c_loc(x_ptr), iexcon_ptr, (/1/))
908 x_ptr => x(iqnstr) ; call c_f_pointer (c_loc(x_ptr), iqnstr_ptr, (/1/))
909
910 call assign (nshl, ndist, nused, irfcon_ptr, iexcon_ptr, iqnstr_ptr, nx, nftw)
911
912 if (ndist == 0) then
913 nstate = 0
914 return
915 end if
916
917 nav = nx - nused
918 nnext = next + nused
919 if (nndel == 0) call print1 (iqnstr_ptr)
920 end if
921
922 if (confpf >= 10 .and. nndel == 0) call print2 (0)
923
924 ! convert some real pointers to integer pointers for use in WFN
925 x_ptr => x(inodi) ; call c_f_pointer (c_loc(x_ptr), nodi_ptr, (/1/))
926 x_ptr => x(indi) ; call c_f_pointer (c_loc(x_ptr), ndi_ptr, (/1/))
927 x_ptr => x(ndel) ; call c_f_pointer (c_loc(x_ptr), ndel_ptr, (/1/))
928 x_ptr => x(iqnstr) ; call c_f_pointer (c_loc(x_ptr), iqnstr_ptr, (/1/))
929 x_ptr => x(nnext) ; call c_f_pointer (c_loc(x_ptr), nnext_ptr, (/1/))
930
931 call wfn (nncsf, nadel, iidist, iidis1, nodi_ptr, ndi_ptr, x(icdi), ndel_ptr, iqnstr_ptr, x(nnext), nnext_ptr, nav)
932
933 if (nndel > 0 .and. nadel > nndel .and. iidis1 == 0) return
934
935 if (nndel /= 0 .and. iidist /= 0) then
936 if (nstate == 1) call print1 (iqnstr_ptr)
937 if (confpf >= 10 .and. iidis1 /= 0) call print2 (iidis1)
938 end if
939
940 ! ascend in coupling tree
941
942 500 do while (kc > 1)
943 kc = kc - 1
944 kcs = kcs - 1
945 qnshl(1,kcs) = qnshl(1,kcs) - 2
946 if (qnshl(1,kcs) >= spnmin(kc)) go to 300
947 end do
948
949 if (nndel /= 0 .and. iidis1 /= 0) call print3 (iidis1, 1)
950 if (nndel == 0) call print3 (idumm, 2)
951
952 end subroutine cplem
953
954
955 !> \brief Distribute electrons to spin-orbitals.
956 !>
957 !> This subroutine will generate all possible configuration state functions subject to setup and constraints set
958 !> by the input namelist. All generated configurations are normally stored in the global arrays \c nodi,
959 !> \c ndi and \c cdi, but if the storage is about to overflow (which is very likely for
960 !> large molecules), it is flushed to the output file sooner (see \ref congen_distribution::wfn), providing space for
961 !> further configurations.
962 !>
963 !> \param nelecp Number of electrons per set of orbitals.
964 !> \param nshlp Number of shells in a set of orbitals.
965 !> \param shlmx Max occupation of a shell or a pseudoshell.
966 !> \param occshl On return, occupation of shells / pseudoshells.
967 !> \param nslsv Save area for index, per shell.
968 !> \param kdssv ?
969 !> \param loopf Pointer for zero occup or pseudoshell, per shell.
970 !> \param ksss Save area for index, per shell.
971 !> \param pqn Shell data (triplet per shell): 1st == 0 (pseudoshell) or /= 0 (real shell index),
972 !> 2nd start index for pseudoshell, 3rd end index for pseudoshell.
973 !> \param occst ?
974 !> \param shlmx1 ?
975 !> \param pqnst ?
976 !> \param mshl Quantum number of a shell or of a pseudoshell, per shell.
977 !> \param mshlst ?
978 !> \param gushl Gerade/ungerade flag, per shell.
979 !> \param gushst ?
980 !> \param cup Coupling scheme.
981 !> \param cupst Expanded shells with zeroes deleted.
982 !> \param ndprod Number of sets of orbitals.
983 !> \param symtyp Symmetry group flag (< 2 diatomic, >= 2 abelian).
984 !> \param confpf ?
985 !> \param sshl ?
986 !> \param sshlst ?
987 !> \param ncsf ?
988 !> \param nncsf ?
989 !> \param nadel ?
990 !> \param nndel ?
991 !> \param x Real workspace.
992 !> \param nftw Unit for text output of the program.
993 !>
994 subroutine distrb (nelecp, nshlp, shlmx, occshl, nslsv, kdssv, loopf, ksss, pqn, occst, shlmx1, pqnst, mshl, mshlst, &
995 gushl, gushst, cup, cupst, ndprod, symtyp, confpf, sshl, sshlst, ncsf, nncsf, nadel, nndel, x, nftw, &
996 noex)
997
998 use precisn, only : wp
999 use congen_data, only : ntyp, ndist, nstate, ncsft => ncsf, confpt => confpf, nshrun => nshl
1000
1001 integer :: confpf, nadel, ncsf, ndprod, nftw, nncsf, nndel, symtyp
1002 integer, dimension(3,*) :: cup, cupst, pqn, pqnst
1003 integer, dimension(*) :: gushl, gushst, kdssv, ksss, loopf, mshl, mshlst, nelecp, nshlp, nslsv, occshl, &
1004 occst, shlmx, shlmx1, sshl, sshlst
1005 real(kind=wp), dimension(*) :: x
1006 intent (in) confpf, cup, gushl, mshl, ncsf, ndprod, nelecp, nndel, nshlp, pqn, shlmx, shlmx1, sshl, symtyp
1007 intent (out) gushst, mshlst, pqnst, sshlst
1008 intent (inout) cupst, kdssv, ksss, loopf, nslsv, occshl, occst
1009
1010 integer :: i, ibias, ic1, id, idd, it, it1, it2, ita, j, kds, kdsb, kdst, kprod, krun, ksi, kss, nadd, nadd2, ncrun, &
1011 nela, neleft, ninitx, nshlw, nslots, noex
1012
1013 confpt = confpf ! Print flag.
1014 ntyp = 0 ! Prototype number ...?
1015 nstate = 0 ! Global number of couplings computed (later) in "cplea" or "cplem".
1016 ndist = 0 ! Number of distributions generated from set of shells (set in "cplea"/"cplem" via "assign").
1017 ibias = 0 ! Offset in the list of shells.
1018 kprod = 0 ! Set of orbitals currently being processed.
1019 ninitx = sum(nshlp(1:ndprod)) ! Total number of shells in all sets of orbitals.
1020
1021 orbital_set_loop: do
1022
1023 kprod = kprod + 1 ! Move on to the next set of orbitals.
1024 kds = 0 ! Shell index within the set of orbitals.
1025 neleft = nelecp(kprod) ! (Remaining) number of movable electrons in the current set of orbitals.
1026 nshlw = nshlp(kprod) ! Number of shells in current set of orbitals.
1027 nslots = sum(shlmx(ibias + 1 : ibias + nshlw)) ! Total number of spin-orbitals in the current set of orbitals.
1028 occshl(ibias + 1 : ibias + nshlw) = 0 ! All shells in this set are un-occupied at the beginning.
1029
1030 shell_loop: do
1031
1032 kds = kds + 1 ! Move on to the next shell (orbital) within the current set of orbitals.
1033 kdsb = kds + ibias ! Global index of the current shell.
1034 occshl(kdsb) = min(neleft, shlmx(kdsb)) ! Put all available electrons in that shell, if possible.
1035
1036 neleft = neleft - occshl(kdsb) ! Electrons that still did not fit into any shell.
1037 nslots = nslots - shlmx(kdsb) ! Number of free spin-orbitals for redistribution of remaining electrons.
1038
1039 if (neleft /= 0) then
1040 if (nslots >= neleft) cycle shell_loop ! Enough room to redistribute the remaining electrons.
1041 do while (nslots <= neleft .and. kds /= 0) ! Repeat until there is either enough space, or all shells are empty.
1042 neleft = neleft + occshl(kdsb) ! Get back electrons from the last populated shell.
1043 occshl(kdsb) = 0 ! Its occupancy is now zero.
1044 nslots = nslots + shlmx(kdsb) ! The number of slots accordingly rises.
1045 kds = kds - 1 ! Return before the last shell.
1046 kdsb = kdsb - 1 ! Return before the last shell (global index).
1047 end do
1048 go to 130
1049 else if (kprod /= ndprod) then ! Not all sets of orbitals processed?
1050 nslsv(kprod) = nslots ! Store the number of free spin-orbitals (of the current set of orbitals).
1051 kdssv(kprod) = kds ! Store the last processed shell (of the current set of orbitals).
1052 ibias = ibias + nshlp(kprod) ! Move the global shell index offset to the beginning of the next set of orbitals.
1053 cycle orbital_set_loop ! Next set of orbitals, please!
1054 end if
1055
1056 ! No electrons left and all sets of orbitals processed here...
1057
1058 where (occshl(1:ninitx) == 0 .or. pqn(1,1:ninitx) /= 0)
1059 loopf(1:ninitx) = 1 ! Un-occupied or movable-to real shells.
1060 elsewhere
1061 loopf(1:ninitx) = 0 ! Occupied and non-movable shells.
1062 end where
1063
1064 ksss(1:ninitx) = 0 !
1065 kdst = 0 !
1066 ksi = 0 ! Orbital set index.
1067
1068 ! ----------------------------------------------------------------------------------------------------------
1069 ! expansion of pseudoshells and delete of zeroes follow the "st" arrays but for cupst are set up
1070
1071 ksi_loop: do
1072
1073 if (ksi == ninitx) then
1074
1075 ! expand and compress the coupling scheme
1076 nshrun = ninitx
1077 ncrun = nshrun - 1
1078 cupst(1:3,1:ncrun) = cup(1:3,1:ncrun)
1079 krun = 0
1080
1081 all_shell_loop: do i = 1, ninitx
1082 krun = krun + 1
1083 if (ksss(i) < 1) then
1084 ! delete section
1085 do id = 1, ncrun
1086 if (cupst(1,id) == krun) then ; it1 = cupst(2,id) ; it2 = cupst(3,id) ; exit ; end if
1087 if (cupst(2,id) == krun) then ; it1 = cupst(1,id) ; it2 = cupst(3,id) ; exit ; end if
1088 end do
1089
1090 ncrun = ncrun - 1
1091
1092 do idd = id, ncrun
1093 cupst(1:3,idd) = cupst(1:3,idd+1)
1094 end do
1095
1096 do idd = 1, ncrun
1097 if (cupst(1,idd) == it2) then ; cupst(1,idd) = it1 ; exit ; end if
1098 if (cupst(2,idd) == it2) then ; cupst(2,idd) = it1 ; exit ; end if
1099 end do
1100
1101 do idd = 1, ncrun
1102 do j = 1, 3
1103 if (cupst(j,idd) <= krun) cycle
1104 cupst(j,idd) = cupst(j,idd) - 1
1105 if (cupst(j,idd) >= it2) cupst(j,idd) = cupst(j,idd) - 1
1106 end do
1107 end do
1108
1109 nshrun = nshrun - 1
1110 krun = krun - 1
1111 else if (ksss(i) > 1) then
1112 ! add section
1113 nadd = ksss(i) - 1
1114 nadd2 = nadd + nadd
1115 do idd = 1, ncrun
1116 cupst(3,idd) = cupst(3,idd) + nadd2
1117 do j = 1, 2
1118 if (cupst(j,idd) == krun) then
1119 cupst(j,idd) = nshrun + nadd2
1120 else if (cupst(j,idd) > krun .and. cupst(j,idd) <= nshrun) then
1121 cupst(j,idd) = cupst(j,idd) + nadd
1122 else if (cupst(j,idd) > krun .and. cupst(j,idd) > nshrun) then
1123 cupst(j,idd) = cupst(j,idd) + nadd2
1124 end if
1125 end do
1126 end do
1127 ic1 = krun
1128 nshrun = nshrun + nadd
1129 do idd = 1, nadd
1130 cupst(1,ncrun+idd) = ic1
1131 cupst(2,ncrun+idd) = krun + idd
1132 cupst(3,ncrun+idd) = nshrun + idd
1133 ic1 = nshrun + idd
1134 end do
1135 ncrun = ncrun + nadd
1136 krun = krun + nadd
1137 end if
1138 end do all_shell_loop ! end of add/delete loop
1139
1140 ! now sort the expanded coupling array
1141 do i = 1, ncrun - 1
1142 do while (i /= cupst(3,i) - kdst)
1143 it = cupst(3,i) - kdst
1144 do j = 1, 3
1145 ita = cupst(j,it)
1146 cupst(j,it) = cupst(j,i)
1147 cupst(j,i) = ita
1148 end do
1149 end do
1150 end do
1151
1152 ! generate a new wave function prototpe
1153 ntyp = ntyp + 1
1154 if (symtyp < 2) then
1155 call cplea (nncsf, nadel, x, nftw)
1156 else
1157 call cplem (nncsf, nadel, x, nftw, noex)
1158 end if
1159 if (nstate == 0) ntyp = ntyp - 1
1160
1161 ! end of coupling scheme section
1162 if (loopf(ksi) == 0) then
1163 go to 320
1164 else
1165 kdst = kdst - ksss(ksi)
1166 go to 360
1167 end if
1168 end if
1169
1170 ksi = ksi + 1 ! Next set of orbitals, please.
1171 kss = 1 ! Shell index.
1172 if (occshl(ksi) == 0) cycle ksi_loop ! Do not process un-occupied shells.
1173 if (occshl(ksi) /= 0) continue
1174
1175 nela = occshl(ksi) ! Occupation of current (ksi) shell.
1176 kss = 0 ! ... ???
1177
1178 410 kss_loop: do while (kss <= pqn(3,ksi) - pqn(2,ksi)) ! Loop over prescribed orbitals in the current set.
1179 kss = kss + 1 ! Next orbital within ksi, please.
1180 kdst = kdst + 1 ! ... ???
1181 occst(kdst) = min(nela, shlmx1(sshl(ksi))) !
1182 if (kss /= 1 .and. kdst > 1) occst(kdst) = min(nela, occst(kdst-1))
1183 pqnst(1:3,kdst) = pqn(1:3,ksi)
1184 mshlst(kdst) = mshl(ksi)
1185 gushst(kdst) = gushl(ksi)
1186 sshlst(kdst) = sshl(ksi)
1187 nela = nela - occst(kdst)
1188 if (nela /= 0) then
1189 cycle kss_loop
1190 else
1191 ksss(ksi) = kss
1192 cycle ksi_loop
1193 end if
1194 end do kss_loop
1195
1196 ! second part of the expansion of pseudoshells and deletion of zeros
1197 320 unkss_loop: do
1198 occst(kdst) = occst(kdst) - 1
1199 nela = nela + 1
1200 if (occst(kdst) /= 0) go to 410
1201 kss = kss - 1
1202 kdst = kdst - 1
1203 if (kss >= 1) then
1204 cycle unkss_loop
1205 else
1206 nela = 0
1207 go to 360
1208 end if
1209 end do unkss_loop
1210
1211 360 unksi_loop: do
1212 ksi = ksi - 1
1213 if (ksi <= 0) exit ksi_loop ! sic
1214 if (loopf(ksi) == 0) then
1215 kss = ksss(ksi)
1216 go to 320
1217 end if
1218 kdst = kdst - ksss(ksi)
1219 end do unksi_loop
1220
1221 ! end of pseudoshell expansion loops
1222
1223 end do ksi_loop
1224
1225 ! ----------------------------------------------------------------------------------------------------------
1226
1227 ! last part of outer loop follows
1228
1229 120 if (occshl(kdsb) == 0) go to 140
1230 occshl(kdsb) = occshl(kdsb) - 1 ! Remove one electron from the current shell.
1231 neleft = 1 ! Remember that we have it
1232 if (nndel /= 0 .and. nadel > nndel) exit orbital_set_loop
1233 if (kds < nshlp(kprod)) cycle shell_loop ! If we can push it into the NEXT shell, do it.
1234 neleft = neleft + occshl(kdsb) ! Otherwise, get back all electrons from the last shell.
1235 occshl(kdsb) = 0 ! Clear its population.
1236 nslots = nslots + shlmx(kdsb) ! Update available number of slots.
1237 kds = kds - 1 ! Move back one shell.
1238 kdsb = kdsb - 1 ! Global index, too.
1239
1240 130 unkds_loop: do while (kds /= 0) ! Reset preceding shells one by one...
1241 if (occshl(kdsb) /= 0) then ! If shell occupied:
1242 occshl(kdsb) = occshl(kdsb) - 1 ! - remove one electron
1243 neleft = neleft + 1 ! - put among the unused electrons
1244 if (nndel /= 0 .and. nadel > nndel) exit orbital_set_loop
1245 cycle shell_loop ! - try to add the electrons to shells
1246 end if
1247 nslots = nslots + shlmx(kdsb) ! If not occupied, add these shells to list of known slots.
1248 kds = kds - 1 ! Recede by one more shell...
1249 kdsb = kdsb - 1 ! Global index too.
1250 end do unkds_loop
1251
1252 140 kprod = kprod - 1 ! Recede by one group of orbitals.
1253 if (kprod == 0) exit orbital_set_loop ! If at beginning, exit subroutine.
1254 ibias = ibias - nshlp(kprod) ! Decrease global shell index offset.
1255 nslots = nslsv(kprod) ! Get (previously stored) number of free slots.
1256 kds = kdssv(kprod) ! Get (previously stored) index of the last shell in that set.
1257 kdsb = kds + ibias ! Get its global index.
1258 go to 120
1259
1260 end do shell_loop
1261
1262 end do orbital_set_loop
1263
1264 ncsft = ncsf
1265
1266 end subroutine distrb
1267
1268
1269 !> \brief Copy determinant data.
1270 !>
1271 !> Whatever was this function originally designed for, it just straghtforwardly copies arrays of data.
1272 !>
1273 !> \param iqn Source data for \c jqn.
1274 !> \param ni Second rank of \c iqn.
1275 !> \param jqn Destination for data from \c iqn.
1276 !> \param nj Second rank of \c jqn.
1277 !> \param ci Source data for \c cj.
1278 !> \param cj Destination for data from \c ci.
1279 !> \param nt Minimal dimension of of \c ci and \c cj and of the third rank of \c iqn and \c jqn.
1280 !>
1281 subroutine packdet (iqn, ni, jqn, nj, ci, cj, nt)
1282
1283 use precisn, only : wp
1284
1285 integer :: ni, nj, nt
1286 real(kind=wp), dimension(*) :: ci, cj
1287 integer, dimension(2,ni,*) :: iqn
1288 integer, dimension(2,nj,*) :: jqn
1289 intent (in) ci, iqn, ni, nj, nt
1290 intent (out) cj, jqn
1291
1292 jqn(1:2,1:nj,1:nt) = iqn(1:2,1:nj,1:nt)
1293 cj(1:nt) = ci(1:nt)
1294
1295 end subroutine packdet
1296
1297
1298 !> \brief ?
1299 !>
1300 subroutine print1 (qnstor)
1301
1302 use global_utils, only : getin
1303 use congen_data , only : confpf, ndist, ne, nshl, ntyp, occshl, pqnr => pqnst, mshl, gushl, qnshl, &
1304 nnelcg => nnlecg, symtyp, nftw, blnk43, header, lp, star, nitem
1305 use congen_pagecontrol, only : addl, newpg, space, taddl, taddl1
1306
1307 integer, dimension(*), intent(in) :: qnstor
1308
1309
1310 integer :: i, ii, imax, j, jj, k, kf, ki, klabel, lt, lta, nlex, nstar
1311 character(len=4), dimension(3,4) :: label = reshape((/ ' NTY', 'P=XX', 'X ', &
1312 ' NDI', 'ST=Y', 'YYY ', &
1313 ' NST', 'ATE=', 'ZZZ ', &
1314 ' ', ' ', ' ' /) , (/ 3, 4/))
1315 character(len=16), parameter :: frmt = '(3A4,A8,I8,8I12)'
1316 character(len=1), dimension(12,4) :: labell
1317
1318 equivalence(label(1,1), labell(1,1))
1319
1320 ! print type(confpf.ge.1) and distribution data (confpf.ge.20)
1321 ne = nnelcg
1322 if (confpf < 1) return
1323 nlex = 0
1324 if (symtyp == 1) nlex = 1
1325 nstar = 20 + min(nitem, nshl) * 12
1326
1327 klabel = 1
1328 call getin (ntyp, 3, labell(7,1), 0)
1329
1330 if (confpf <= 1 .and. ntyp /= 1) then
1331 call taddl1 (5 + (6 + nlex) * ((nshl + nitem - 1) / nitem), lta)
1332 if (lta /= 0) then
1333 call space (2)
1334 else
1335 call newpg
1336 end if
1337 else
1338 call newpg
1339 end if
1340
1341 do i = 1, nshl, nitem
1342 imax = min(i + nitem - 1, nshl)
1343 call addl (5 + nlex)
1344 write(nftw,frmt) (label(j,klabel), j=1,3), header(1), (j, j=i,imax)
1345 write(nftw,frmt) blnk43, header(2), (occshl(j), j=i,imax)
1346 write(nftw,frmt) blnk43, header(3), (mshl(j), j=i,imax)
1347 if (symtyp == 1) write(nftw,frmt) blnk43, header(4), (gushl(j), j=i,imax)
1348
1349 ! print occshl and sym data
1350
1351 klabel = 4
1352 write(nftw,'(3A4,A8,9(A3,2(I3,","),I3,")"))') blnk43, header(5), (lp, (pqnr(jj,j), jj=1,3), j=i,imax)
1353 if (symtyp < 2) then
1354 write(nftw,'(3A4,A8,9(A3,2(I2,","),I2,")"))') blnk43, header(6), (lp, (qnshl(jj,j), jj=1,3), j=i,imax)
1355 else
1356 write(nftw,frmt) blnk43, header(6), (qnshl(1,j), j=i,imax)
1357 end if
1358 call space (1)
1359 end do
1360
1361 ! bypass printing of distributions if not required
1362
1363 if (confpf >= 20) then
1364
1365 ! print row of star separator
1366
1367 call taddl (1, lt)
1368 if (lt > 0) write(nftw,'(" ",132A1)') (star, j=1,nstar)
1369 call space (1)
1370 kf = 0
1371 do ii = 1, ndist
1372 klabel = 2
1373 call getin (ii, 4, labell(8,2), 0)
1374 do i = 1, nshl, nitem
1375 imax = min(i + nitem - 1, nshl)
1376 call addl (4 + nlex)
1377 write(nftw,frmt) (label(j,klabel), j=1,3), header(1), (j, j=i,imax)
1378 write(nftw,frmt) blnk43, header(2), (occshl(j), j=i,imax)
1379 write(nftw,frmt) blnk43, header(3), (mshl(j), j=i,imax)
1380 if (symtyp == 1) write(nftw,frmt) blnk43, header(4), (gushl(j), j=i,imax)
1381 klabel = 4
1382 ki = kf + 1
1383 kf = ki + (imax - i)
1384 write(nftw,frmt) blnk43, header(5), (qnstor(k), k=ki,kf)
1385 call space (1)
1386 end do
1387 end do
1388
1389 end if
1390
1391 ! merge confpf >= 1 and confpf >= 20 paths
1392
1393 call addl (1)
1394 write(nftw,'(" ",19("*"),5X,"TOTAL NUMBER OF DISTRIBUTIONS FOR NTYP =",I3," IS",I5)') ntyp, ndist
1395 if (confpf < 10) return
1396
1397 ! prepare for state printing
1398
1399 call space (1)
1400 call taddl (1, lt)
1401 if (lt > 0) write(nftw,'(" ",132A1)') (star, j=1,nstar)
1402 call space (1)
1403
1404 end subroutine print1
1405
1406
1407 !> \brief ?
1408 !>
1409 subroutine print2 (iidis1)
1410
1411 use global_utils, only : getin
1412 use congen_data, only : confpf, ncsf, ndist, nshl, nstate, cup, qnshl, symtyp, nftw, blnk43, header, lp, nitem
1413 use congen_pagecontrol, only : addl, space
1414
1415 integer :: iidis1
1416 intent (in) iidis1
1417
1418 integer :: i, imax, j, jj, kf, ki, klabel, ncsff, ncsfi, nshlm1
1419 character(len=4), dimension(3,4) :: label = reshape((/' NTY', 'P=XX', 'X ', &
1420 ' NDI', 'ST=Y', 'YYY ', &
1421 ' NST', 'ATE=', 'ZZZ ', &
1422 ' ', ' ', ' '/) , (/ 3, 4/))
1423 character(len=1), dimension(12,4) :: labell
1424
1425 equivalence(label(1,1), labell(1,1))
1426
1427 ! print state data
1428
1429 nshlm1 = nshl - 1
1430 if (confpf < 10) return
1431 call getin (nstate, 3, labell(9,3), 0)
1432 if (nshlm1 <= 0) then
1433 call addl (1)
1434 write(nftw,'(3A4,A8,I8,8I12)') (label(j,3), j=1,3)
1435 return
1436 end if
1437 kf = nshl
1438 klabel = 3
1439 do i = 1, nshlm1, nitem
1440 imax = min(i + nitem - 1, nshlm1)
1441 call addl (2) ! JMC argument should probably be 3 as there are 3 writes below???
1442 write(nftw,'(3A4,A8,9(A3,2(I2,","),I2,")"))') (label(j,klabel), j=1,3), header(7), (lp, (cup(jj,j), jj=1,3), j=i,imax)
1443 klabel = 4
1444 ki = kf + 1
1445 kf = ki + (imax - i)
1446 write(nftw,'(3A4,A8,9(A3,2(I2,","),I2,")"))') blnk43, header(6), (lp, (qnshl(jj,j), jj=1,3), j=ki,kf)
1447 if (symtyp < 2) then
1448 write(nftw,'(3A4,A8,9(A3,2(I2,","),I2,")"))') blnk43, header(6), (lp, (qnshl(jj,j), jj=1,3), j=ki,kf)
1449 else
1450 write(nftw,'(3A4,A8,I8,8I12)') blnk43, header(6), (qnshl(1,j), j=ki,kf)
1451 end if
1452 call space (1)
1453 end do
1454
1455 ! print CSF numbers for this state
1456
1457 if (iidis1 /= 0) then
1458 ncsfi = ncsf - iidis1 + 1
1459 ncsff = ncsf
1460 else
1461 ncsfi = ncsf + 1
1462 ncsff = ncsf + ndist
1463 end if
1464
1465 call addl (1)
1466 write(nftw,'(" ",19("*"),5X,"CSF NUMBERS",I6," TO",I6," GENERATED FOR NSTATE=",I3)') ncsfi, ncsff, nstate
1467 call space (1)
1468
1469 end subroutine print2
1470
1471
1472 !> \brief ?
1473 !>
1474 subroutine print3 (iidis1, i13)
1475
1476 use congen_data, only : confpf, ncsf, ndist, nshl, nstate, ntyp, nftw, star, nitem
1477 use congen_pagecontrol, only : addl, space, taddl
1478
1479 integer, intent(in) :: i13, iidis1
1480
1481 integer :: i, j, lt, ncsfi, nstar, nstot
1482
1483 ! print summary data
1484 if (i13 /= 2) then
1485 nstot = iidis1
1486 else
1487 if (confpf < 1) return
1488 nstot = nstate * ndist
1489 end if
1490 ncsfi = ncsf + 1 - nstot
1491 if (confpf >= 40) call space (1)
1492 call addl (2)
1493
1494 write(nftw,'(" ",19("*"),5X,"TOTAL NUMBER OF STATES FOR NTYP=",I3," IS",I4)') ntyp, nstate
1495 write(nftw,'(20X,"CSF NUMBERS",I10," TO",I10," (",I9," CSFS) GENERATED")') ncsfi, ncsf, nstot
1496
1497 nstar = 20 + min(nitem, nshl) * 12 ! JMC see PRINT1 to understand where this hardwiring has come from
1498 do i = 1, 2
1499 call taddl (1, lt)
1500 if (lt > 0) write(nftw,'(" ",132A1)') (star, j=1,nstar)
1501 end do
1502
1503 end subroutine print3
1504
1505
1506 !> \brief ?
1507 !>
1508 subroutine print4 (nd,x,ix)
1509
1510 use precisn, only : wp
1511 use congen_data, only : ne, nftw
1512 use congen_pagecontrol, only : addl, space
1513
1514 integer :: nd
1515 integer, dimension(*) :: ix
1516 real(kind=wp), dimension(*) :: x
1517 intent (in) ix, nd, x
1518
1519 integer :: ie, ind, ip, ipi, it, j, nl, nitem
1520
1521 ! print determinant and spin-orbital data
1522
1523 ipi = 0
1524 it = 1
1525 nitem = 20
1526 nl = (ne + nitem - 1) / nitem
1527
1528 do ind = 1, nd
1529 call addl (nl + it)
1530 if (it /= 0) then
1531 write(nftw,'(30X,"NUMBER OF DET IS",I3)') nd
1532 it = 0
1533 end if
1534 do ie = 1, ne, nitem
1535 ip = ipi + 1
1536 ipi = ip + min(nitem - 1, ne - ie)
1537 if (ie /= 1) then
1538 write(nftw,'(45X,20I4)') (ix(j), j=ip,ipi)
1539 else
1540 write(nftw,'(I25,5X,E15.8,20I4)') ind, x(ind), (ix(j), j=ip,ipi)
1541 end if
1542 end do
1543 end do
1544
1545 call space (1)
1546
1547 end subroutine print4
1548
1549
1550 !> \brief ?
1551 !>
1552 subroutine print5 (nd, ndi)
1553
1554 use congen_data, only : ncsf, nftw
1555 use congen_pagecontrol, only : addl
1556
1557 integer :: nd
1558 integer, dimension(*) :: ndi
1559 intent (in) nd, ndi
1560
1561 integer :: ind, ip, ipi, j, nrep
1562
1563 ip = 1
1564 do ind = 1, nd
1565 nrep = ndi(ip)
1566 if (nrep == 0) then
1567 call addl (1)
1568 write(nftw, '(35X,2I5)') ind, nrep
1569 if (ind == 1) then
1570 write(nftw,'("+",29X,I5)') ncsf
1571 end if
1572 else
1573 call addl (2)
1574 ipi = ip + 1
1575 ip = ip + nrep
1576
1577 if (ind /= 1) then
1578 write(nftw,'(35X,2I5,21I4/(45X,21I4))') ind, nrep, (ndi(j), j=ipi,ip)
1579 else
1580 write(nftw,'(30X,3I5,21I4/(45X,21I4))') ncsf, ind, nrep, (ndi(j), j=ipi,ip)
1581 end if
1582
1583 ipi = ip + 1
1584 ip = ip + nrep
1585
1586 write(nftw,'(45X,21I4)') (ndi(j), j = ipi,ip)
1587 end if
1588 ip = ip + 1
1589 end do
1590
1591 end subroutine print5
1592
1593
1594 !> \brief ?
1595 !>
1596 !> \param ns Actually, this is just reference to the global \ref congen_data::nshl (number of shells in the set of orbitals
1597 !> currently being processed.)
1598 !> \param x Real workspace (or the currently free part of it).
1599 !> \param last Number of available elements in the workspace.
1600 !> \param nd ?
1601 !> \param confpf Print level flag.
1602 !>
1603 subroutine state (ns, x, last, nd, confpf)
1604
1605 use iso_c_binding, only : c_loc, c_f_pointer
1606 use precisn, only : wp
1607 use congen_data, only : lratio, cup, iqn => qnshl, ne => nnlecg, isz => nisz
1608
1609 integer :: confpf, last, nd, ns
1610 real(kind=wp), dimension(*), target :: x
1611 intent (in) confpf, last
1612 intent (inout) nd, x
1613
1614 integer :: ic, id, intpfg, iqns, jqns, lc, ld, ldp, ldq, nam, ndmx, nt, ntmx
1615
1616 real(wp), pointer :: x_ptr
1617 integer, pointer, dimension(:) :: iqns_ptr, jqns_ptr, ld_ptr
1618
1619 x_ptr => x(1)
1620
1621 intpfg = merge(1, 0, confpf > 40)
1622
1623 nd = 0
1624 nam = ns + ns - 1
1625 ntmx = last / (nam + nam + 1)
1626 ic = 1
1627 iqns = ic + ntmx
1628
1629 ! convert real pointer to integer pointer for use in WFCPLE
1630 x_ptr => x(iqns) ; call c_f_pointer (c_loc(x_ptr), iqns_ptr, (/1/))
1631
1632 call wfcple (nam, iqn, isz, cup, iqns_ptr, x(ic), ntmx, nt, intpfg)
1633
1634 if (nt == 0) return
1635
1636 jqns = last - 2 * ns * nt + 1
1637 lc = jqns - nt
1638
1639 ! convert real pointer to integer pointer for use in PACK
1640 x_ptr => x(jqns) ; call c_f_pointer (c_loc(x_ptr), jqns_ptr, (/1/))
1641
1642 call packdet (iqns_ptr, nam, jqns_ptr, ns, x(1), x(lc), nt)
1643
1644 ndmx = (lc - 1) / (ne + 1)
1645 ld = 1 + ndmx
1646
1647 ! convert real pointer to integer pointer for use in GETSO
1648 x_ptr => x(ld) ; call c_f_pointer (c_loc(x_ptr), ld_ptr, (/1/))
1649
1650 call getso (ns, intpfg, nt, jqns_ptr, x(lc), nd, ld_ptr, x(ic), ndmx)
1651
1652 if (nd == 0) return
1653
1654 ! Finalize
1655
1656 ldp = nd + 1
1657 ldq = nd + (ne * nd + lratio - 1) / lratio
1658 do id = ldp, ldq
1659 x(id) = x(ld)
1660 ld = ld + 1
1661 end do
1662
1663 end subroutine state
1664
1665
1666 !> \brief ?
1667 !>
1668 subroutine wfcple (nam, iqn, isz, icup, iqns, c, last, lc2, intpfg)
1669
1670 use precisn, only : wp
1671 use consts, only : one => xone
1672 use congen_data, only : root2, nftw, lg
1673
1674 integer :: intpfg, isz, last, lc2, nam
1675 real(kind=wp), dimension(*) :: c
1676 integer, dimension(3,lg) :: icup
1677 integer, dimension(3,*) :: iqn
1678 integer, dimension(2,nam,*) :: iqns
1679 intent (in) icup, isz, last, nam
1680 intent (inout) c, iqns, lc2
1681
1682 integer :: i, iam, ic, init, j, l, lc, lc1, lc3, m, mp, ms, n, n1, n2, n3, nc, niam, niam1
1683 logical, dimension(nam) :: ind ! JMC changing the dimension from 150
1684 integer, dimension(2,200) :: iszt ! JMC not sure how to dimension this...
1685 real(kind=wp) :: sign
1686
1687 iqns(1,nam,1) = isz
1688 iqns(2,nam,1) = iqn(2,nam)
1689 c(1) = one
1690 lc2 = 1
1691 if (nam == 1) return
1692 do i = 1, nam
1693 ind(i) = .true.
1694 end do
1695 niam = (nam + 1) / 2
1696 100 n3 = nam
1697 110 if (ind(n3)) then
1698 if (n3 <= niam) go to 299
1699 ind(n3) = .false.
1700 init = last - lc2 + 1
1701 l = last
1702 lc = lc2
1703 do while (lc > 0)
1704 iqns(1:2,1:nam,l) = iqns(1:2,1:nam,lc)
1705 c(l) = c(lc)
1706 lc = lc - 1
1707 l = l - 1
1708 end do
1709 n = 1
1710 300 if (icup(3,n) == n3) then
1711 n1 = icup(1,n)
1712 n2 = icup(2,n)
1713 if (n1 > n3 .or. n2 > n3) go to 299
1714 if (.not.ind(n1) .or. .not.ind(n2)) go to 299
1715 if (n1 <= niam) ind(n1) = .false.
1716 if (n2 <= niam) ind(n2) = .false.
1717 lc2 = 0
1718 do l = init, last
1719 iqns(2,n1,l) = iqn(2,n1)
1720 iqns(2,n2,l) = iqn(2,n2)
1721 m = iqns(2,n3,l)
1722 if (abs(m) /= iqn(2,n1) + iqn(2,n2)) then
1723 mp = iqn(2,n1) - iqn(2,n2)
1724 if (abs(m) /= abs(mp)) go to 511
1725 n = n2
1726 if (m /= mp) n = n1
1727 iqns(2,n,l) = -iqns(2,n,l)
1728 else if (iqns(2,n3,l) < 0) then
1729 iqns(2,n1,l) = -iqns(2,n1,l)
1730 iqns(2,n2,l) = -iqns(2,n2,l)
1731 end if
1732 lc1 = lc2 + 1
1733 ms = iqns(1,n3,l)
1734 call cgcoef (iqn(1,n1), iqn(1,n2), iqn(1,n3), ms, nc, iszt, c(lc1), intpfg)
1735 if (nc > 0) then
1736 lc2 = lc2 + nc
1737 if (lc2 >= l) go to 899
1738 ic = 1
1739 do lc = lc1, lc2
1740 iqns(1:2,1:nam,lc) = iqns(1:2,1:nam,l)
1741 iqns(1,n1,lc) = iszt(1,ic)
1742 iqns(1,n2,lc) = iszt(2,ic)
1743 ic = ic + 1
1744 c(lc) = c(lc) * c(l)
1745 end do
1746 if (iqns(2,n3,l) < 0) then
1747 if (iqn(3,n1) >= 0 .and. iqn(3,n2) >= 0) cycle
1748 c(lc1:lc2) = -c(lc1:lc2)
1749 else if (iqns(2,n3,l) == 0) then
1750 if (iqn(2,n1) /= 0) then
1751 if (lc2 + nc >= l) go to 899
1752 lc3 = lc2
1753 sign = one
1754 if (iqn(3,n3) < 0) sign = -one
1755 do lc = lc1, lc3
1756 lc2 = lc2 + 1
1757 iqns(1:2,1:nam,lc2) = iqns(1:2,1:nam,lc)
1758 iqns(2,n1,lc2) = -iqns(2,n1,lc2)
1759 iqns(2,n2,lc2) = -iqns(2,n2,lc2)
1760 c(lc) = c(lc) * root2
1761 c(lc2) = c(lc) * sign
1762 end do
1763 cycle
1764 end if
1765 if (iqn(3,n1) * iqn(3,n2) /= iqn(3,n3)) go to 511
1766 end if
1767 cycle
1768 end if
1769 go to 511
1770 end do
1771 go to 100
1772 end if
1773 n = n + 1
1774 if (n < niam) go to 300
1775 go to 299
1776 end if
1777 n3 = n3 - 1
1778 if (n3 > 0) go to 110
1779 return
1780
1781 299 niam1 = niam - 1
1782 write(nftw,'("0ERROR IN COUPLING TREE"//(3I5))') ((icup(i,j), i=1,3), j=1,niam1)
1783 lc2 = 0
1784 return
1785
1786 511 write(nftw,'("0COUPLING IMPOSSIBLE MULTIPLICITIES FOLLOW"//(3I5))') (iqn(i,n1), iqn(i,n2), iqn(i,n3), i=1,3)
1787 lc2 = 0
1788 return
1789
1790 899 write(nftw,'("0STORAGE OVERFLOW IN VECTOR COUPLING")')
1791 lc2 = 0
1792
1793 end subroutine wfcple
1794
1795
1796 !> \brief TODO...
1797 !>
1798 !> One of interesting features of this subroutine is that it can detect if the storage for determinants is about to
1799 !> overflow. In such a case it will automatically flush the data to the output file (not waiting for the call to
1800 !> \ref congen_driver::csfout), providing space for further calculation.
1801 !>
1802 !> \param nncsf Number of CSFs (total).
1803 !> \param nadel Number of processed configurations.
1804 !> \param iidist ?
1805 !> \param iidis3 ?
1806 !> \param nodi Number of determinants per CSF (1 ... \c noi).
1807 !> \param ndi Compressed determinants (1 ... \c nid).
1808 !> \param cdi Determinant coefficients (1 ... \c nid).
1809 !> \param ndel Configurations read from input.
1810 !> \param pqnshl ?
1811 !> \param x Real workspace, far enough after the chunk corresponding to \c pqnshl.
1812 !> \param ix Actually, the same address as \c x, so the same workspace (just cast to integer).
1813 !> \param nx Available elements in the real workspace.
1814 !>
1815 subroutine wfn (nncsf, nadel, iidist, iidis3, nodi, ndi, cdi, ndel, pqnshl, x, ix, nx)
1816
1817 ! noi # of states (nodi)
1818 ! nid # of determinants (cdi)
1819 ! ni # of replacements (ndi)
1820
1821 use precisn, only : wp
1822 use congen_data, only : lratio, iidis2, lcdi, lndi, ni, nid, noi, exdet, exref, norep, nonew, ntso, &
1823 noimx => nodimx, nidmx => cdimx, jmx => ndimx, megul, nsoi, ncall, nftw, &
1824 confpf, ncsf, ndist, nshl, occshl, sshl => sshlst, shlmx1, nndel, nelec => nnlecg, ndimx
1825
1826 integer :: iidis3, iidist, nadel, nncsf, nx
1827 real(kind=wp), dimension(*) :: cdi, x
1828 integer, dimension(*) :: ix, ndel, ndi, nodi, pqnshl
1829 intent (in) ndel, pqnshl
1830 intent (out) iidis3
1831 intent (inout) cdi, iidist, nadel, ndi, nncsf, nodi
1832
1833 integer :: det, i, ia, ib, id, idist, iidis1, irep, j, k, k1, kshl, kshlst, lf, li, nd, nii
1834
1835 save iidis1 ! JMC adding this in consultation with Jonathan Tennyson because the variable was found to be used but not set in tests...
1836 ! The other variables initialized in the 1st call to this routine are (former) common block variables now in congen_data.
1837
1838 ! reset counters for a new wave function group ("ncall = 1" is set in CSFGEN before calling DISTRB)
1839 if (ncall == 1) then
1840 lcdi = 0
1841 lndi = 0
1842 noi = 0
1843 nid = 0
1844 iidis1 = 0
1845 iidis2 = 0
1846 ni = 0
1847 ncall = 0
1848 end if
1849
1850 ! here most work is done
1851 call state (nshl, x, nx, nd, confpf)
1852
1853 if (nd <= 0) then
1854 write(nftw,'("1","******* ERROR IN WFN ND=0"//)')
1855 stop 70
1856 end if
1857
1858 ! no distributions generated in ASSIGN
1859 if (ndist == 0) return
1860
1861 if (nndel /= 0 .and. nadel > nndel) return
1862
1863 iidist = 0
1864 kshlst = 0
1865
1866 ! for all distributions generated by ASSIGN
1867 do idist = 1, ndist
1868
1869 nncsf = nncsf + 1
1870
1871 if (nndel /= 0 .and. ndel(nadel) /= nncsf) then
1872 kshlst = kshlst + nshl
1873 cycle
1874 end if
1875
1876 ncsf = ncsf + 1
1877 if (nndel /= 0) then
1878 nadel = nadel + 1
1879 iidist = iidist + 1
1880 end if
1881
1882 k = nd * lratio
1883 ia = nd * (2 * nelec + 1) + ni
1884 ib = nd + nid
1885 k1 = noi + 1
1886
1887 if (ia > jmx .or. ib > nidmx .or. k1 > noimx) then
1888 if (nndel == 0 .or. iidis2 /= 0) then
1889 write(megul) noi, (nodi(i), i=1,noi)
1890 write(megul) nid, (cdi(i), i=1,nid)
1891 write(megul) ni, (ndi(i), i=1,ni)
1892 lcdi = lcdi + nid
1893 lndi = lndi + ni
1894 end if
1895 noi = 0
1896 ni = 0
1897 nid = 0
1898 end if
1899
1900 nii = ni + 1
1901 do id = 1, nd
1902 exdet(1:ntso) = 0
1903 kshl = kshlst
1904 lf = 0
1905 do i = 1, nshl
1906 li = lf + 1
1907 lf = lf + occshl(i)
1908 kshl = kshl + 1
1909 do j = li, lf
1910 k = k + 1
1911 det = ix(k) + nsoi(sshl(i)) + (pqnshl(kshl) - 1) * shlmx1(sshl(i))
1912 exdet(det) = 1
1913 end do
1914 end do
1915 ia = 1
1916 ib = 1
1917 do i = 1, ntso
1918 if (exdet(i) < exref(i)) then
1919 norep(ia) = i
1920 ia = ia + 1
1921 else if (exdet(i) > exref(i)) then
1922 nonew(ib) = i
1923 ib = ib + 1
1924 end if
1925 end do
1926 irep = ib - 1
1927 ni = ni + 1
1928
1929 ! ZM check that NI does not overflow NDIMX
1930 ! If NDI was allocatable we could resize it here
1931 if (ni + 2 * irep > ndimx) then
1932 write(nftw,'("******* ERROR IN WFN, NDIMX TOO SMALL",2I10)') ni + 2 * irep, ndimx
1933 stop 81
1934 end if
1935
1936 ndi(ni) = irep
1937 do i = 1, irep
1938 ni = ni + 1
1939 ndi(ni) = norep(i)
1940 ndi(ni + irep) = nonew(i)
1941 end do
1942 ni = ni + irep
1943
1944 end do
1945
1946 if (confpf >= 40) call print5 (nd, ndi(nii))
1947 noi = noi + 1
1948
1949 ! ZM check that NOI does not overflow NODIMX
1950 ! If NODI was allocatable we could resize it here
1951 if (noi > noimx) then
1952 write(nftw,'("******* ERROR IN WFN, NODIMX TOO SMALL",2I10)') noi, noimx
1953 stop 82
1954 end if
1955
1956 nodi(noi) = nd
1957
1958 ! ZM check that NID+ND does not overflow CDIMX
1959 ! If CDI was allocatable we could resize it here
1960 if (nid + nd > nidmx) then
1961 write(nftw,'("******* ERROR IN WFN, CDIMX TOO SMALL",2I10)') nid + nd, nidmx
1962 stop 83
1963 end if
1964
1965 do i = 1, nd
1966 nid = nid + 1
1967 cdi(nid) = x(i)
1968 end do
1969
1970 kshlst = kshl
1971
1972 end do
1973
1974 iidis1 = iidis1 + iidist
1975 iidis2 = iidist + iidis2
1976 if (confpf >= 30 .and. iidist > 0) call print4 (nd, x, ix(nd*lratio+1))
1977 iidis3 = iidis1
1978
1979 end subroutine wfn
1980
1981end module congen_distribution
congen_data
Module containing parameters / data required throughout the CONGEN program.
Definition congen_data.f90:32
congen_data::ndel
integer ndel
Workspace position used for the CSFs read from input.
Definition congen_data.f90:81
congen_data::nobt
integer nobt
Total number of all orbitals considered in the calculation, essentially equal to sum(nob).
Definition congen_data.f90:110
congen_data::root2
real(kind=wp), parameter root2
Note that this is actually inverse square root!
Definition congen_data.f90:51
congen_data::shlmx1
integer, dimension(nu) shlmx1
Maximal occupancy of orbitals per symmetry.
Definition congen_data.f90:103
congen_data::lndi
integer lndi
Total number of integers forming packed determinants.
Definition congen_data.f90:143
congen_data::norep
integer, dimension(lg) norep
Definition congen_data.f90:138
congen_data::test
integer, dimension(jz) test
Determines which of the two types of constraints will be used.
Definition congen_data.f90:117
congen_data::icdi
integer icdi
Position in workspace where the determinant coefficients start (one coeff per one determinant).
Definition congen_data.f90:75
congen_data::gushl
integer, dimension(lg) gushl
Gerade/ungerade per shell.
Definition congen_data.f90:92
congen_data::nst
integer, dimension(nu, lg) nst
Pointer to shell index and order count; only used in congen_distribution::assign.
Definition congen_data.f90:118
congen_data::gutot
integer gutot
Total gerade (= +1) / ungerade (= -1) quantum number.
Definition congen_data.f90:84
congen_data::iqnstr
integer iqnstr
Starting position in workspace for per-shell quantum numbers (e.g. congen_distribution::cplem).
Definition congen_data.f90:79
congen_data::ne
integer ne
Definition congen_data.f90:122
congen_data::indi
integer indi
Position in workspace where packed determinants start (each as size + replacements).
Definition congen_data.f90:77
congen_data::binom
real(kind=wp), dimension((jsmax *jsmax+3 *jsmax+4)/2) binom
Table of precomputed Pascal triangle.
Definition congen_data.f90:106
congen_data::nisz
integer nisz
Total S_z, set by CSFGEN, defaults to first element of qntot.
Definition congen_data.f90:85
congen_data::blnk43
character(len=4), dimension(3), parameter blnk43
Definition congen_data.f90:59
congen_data::qnshlr
integer, dimension(lg) qnshlr
Total quantum numbers for all shells in the set of orbitals being processed.
Definition congen_data.f90:95
congen_data::ntyp
integer ntyp
Prototype number, updated in congen_distribution::distrb.
Definition congen_data.f90:124
congen_data::megul
integer megul
File unit for binary output of generated configurations.
Definition congen_data.f90:128
congen_data::ntso
integer ntso
Total number of spin-orbitals, given the orbitals and their maximal occupacy (from group properties).
Definition congen_data.f90:133
congen_data::nitem
integer, parameter nitem
Number of interesting things printed per line in the PRINT* routines.
Definition congen_data.f90:55
congen_data::exref
integer, dimension(:), allocatable exref
Population of spin-orbitals (0 = not populated, 1 = populated).
Definition congen_data.f90:136
congen_data::mclim
integer, dimension(2, lg) mclim
Definition congen_data.f90:99
congen_data::qntot
integer, dimension(3) qntot
Total quantum numbers.
Definition congen_data.f90:102
congen_data::nrcon
integer, dimension(jz) nrcon
Definition congen_data.f90:116
congen_data::nstate
integer nstate
Number of couplings computed in congen_distribution::cplea or congen_distribution::cplem.
Definition congen_data.f90:127
congen_data::nndel
integer nndel
Number of workspace elements used for reservation of space needed when reading CSFs from input.
Definition congen_data.f90:86
congen_data::ind
integer, dimension(jsmax+2) ind
Pascal triangle row pointers into congen_data::binom.
Definition congen_data.f90:107
congen_data::nshl
integer nshl
Number of shells in the set of orbitals currently being processed.
Definition congen_data.f90:123
congen_data::exdet
integer, dimension(:), allocatable exdet
Definition congen_data.f90:135
congen_data::confpf
integer confpf
Print flag for the amount of print given of configurations and prototypes.
Definition congen_data.f90:120
congen_data::kslim
integer, dimension(2, lg) kslim
Definition congen_data.f90:98
congen_data::nnlecg
integer nnlecg
The total number of electrons which are movable within the current wfn group.
Definition congen_data.f90:87
congen_data::nsoi
integer, dimension(nu) nsoi
Running index of the first spin-orbital within given total symmetry.
Definition congen_data.f90:139
congen_data::ift
integer ift
Used to signal to congen_distribution::assign to initialize some variables.
Definition congen_data.f90:109
congen_data::ndist
integer ndist
Number of distributions generated from set of shells (set in congen_distribution::assign).
Definition congen_data.f90:126
congen_data::nonew
integer, dimension(lg) nonew
Definition congen_data.f90:137
congen_data::nncsf
integer nncsf
Number of CSFs generated by congen_distribution::wfn (total).
Definition congen_data.f90:131
congen_data::lp
character(len=3), parameter lp
Definition congen_data.f90:66
congen_data::pqnst
integer, dimension(3, lg) pqnst
Definition congen_data.f90:91
congen_data::nobi
integer, dimension(nu) nobi
Running index of the first orbital in each symmetry.
Definition congen_data.f90:114
congen_data::lratio
integer lratio
Ratio of real size to integer size. Used to manage workspace data.
Definition congen_data.f90:74
congen_data::spnmin
integer, dimension(lg) spnmin
Definition congen_data.f90:96
congen_data::cdimx
integer cdimx
Maximal number of determinants.
Definition congen_data.f90:125
congen_data::nid
integer nid
Length of the integer array containing all packed determinants.
Definition congen_data.f90:145
congen_data::qnshl
integer, dimension(3, 2 *lg) qnshl
Definition congen_data.f90:100
congen_data::ndimx
integer ndimx
Maximal number of workspace elements usable for packed determinant data.
Definition congen_data.f90:130
congen_data::ncall
integer ncall
Used to signal to "wfn" to initialize some variables.
Definition congen_data.f90:129
congen_data::iexcon
integer iexcon
Position in workspace where EX constraints start.
Definition congen_data.f90:76
congen_data::lcdi
integer lcdi
Total number of determinants (including those already flushed to disk).
Definition congen_data.f90:142
congen_data::irfcon
integer irfcon
Definition congen_data.f90:80
congen_data::nshsym
integer, dimension(nu) nshsym
Definition congen_data.f90:115
congen_data::cup
integer, dimension(3, lg) cup
Coupling scheme data. ???
Definition congen_data.f90:90
congen_data::ncsf
integer ncsf
Total number of CSFs generated.
Definition congen_data.f90:121
congen_data::nftw
integer nftw
Unit number for printing; may be changed via the STATE namelist so not a parameter.
Definition congen_data.f90:57
congen_data::next
integer next
Position in the workspace indicating the first unused element.
Definition congen_data.f90:147
congen_data::thresh1
real(kind=wp), parameter thresh1
Threshold used in congen_distribution::cgcoef.
Definition congen_data.f90:52
congen_data::iidis2
integer iidis2
Definition congen_data.f90:141
congen_data::star
character(len=1), parameter star
Definition congen_data.f90:67
congen_data::noi
integer noi
Number of determinants per CSF.
Definition congen_data.f90:146
congen_data::nx
integer nx
Full or remaining workspace size (depending on context).
Definition congen_data.f90:148
congen_data::header
character(len=8), dimension(7), parameter header
Definition congen_data.f90:61
congen_data::mshl
integer, dimension(lg) mshl
Symmetry number from zero to n-1 (mvalue).
Definition congen_data.f90:93
congen_data::qntar
integer, dimension(3) qntar
Coupling control for current wfn group. ???
Definition congen_data.f90:101
congen_data::nodimx
integer nodimx
Maximal number of workspace elements usable for CSF information.
Definition congen_data.f90:132
congen_data::ntcon
integer ntcon
Definition congen_data.f90:111
congen_data::ni
integer ni
Number of determinants currently in memory.
Definition congen_data.f90:144
congen_data::lg
integer, parameter lg
??? used to dimension arrays below; also in CPLE[A,M].
Definition congen_data.f90:39
congen_data::sshlst
integer, dimension(lg) sshlst
Definition congen_data.f90:97
congen_data::inodi
integer inodi
Position in workspace where CSFs start (as number of determinants per CSF).
Definition congen_data.f90:78
congen_data::symtyp
integer symtyp
Symmetry type, 0 = C_infv, 1 = D_infh, 2 = {D_2h, C_2v, C_s, E}.
Definition congen_data.f90:88
congen_data::occshl
integer, dimension(lg) occshl
Compressed shell occupations, all zeros deleted and pseudo shells expanded.
Definition congen_data.f90:94
congen_distribution
Distribute electrons to available orbitals.
Definition congen_distribution.f90:23
congen_distribution::wfcple
subroutine wfcple(nam, iqn, isz, icup, iqns, c, last, lc2, intpfg)
?
Definition congen_distribution.f90:1669
congen_distribution::distrb
subroutine, public distrb(nelecp, nshlp, shlmx, occshl, nslsv, kdssv, loopf, ksss, pqn, occst, shlmx1, pqnst, mshl, mshlst, gushl, gushst, cup, cupst, ndprod, symtyp, confpf, sshl, sshlst, ncsf, nncsf, nadel, nndel, x, nftw, noex)
Distribute electrons to spin-orbitals.
Definition congen_distribution.f90:997
congen_distribution::print3
subroutine, private print3(iidis1, i13)
?
Definition congen_distribution.f90:1475
congen_distribution::cplem
subroutine, private cplem(nncsf, nadel, x, nftw, noex)
?
Definition congen_distribution.f90:802
congen_distribution::print1
subroutine, private print1(qnstor)
?
Definition congen_distribution.f90:1301
congen_distribution::cplea
subroutine, private cplea(nncsf, nadel, x, nftw)
Loop through (and fill) all allowed couplings for a given electron distribution into shells.
Definition congen_distribution.f90:575
congen_distribution::getsm
subroutine, private getsm(ne, isz, nc, c, iso)
Form ne electrons in a shell coupled to (l,is).
Definition congen_distribution.f90:375
congen_distribution::getso
subroutine, private getso(ns, intpfg, nti, iqns, ci, nd, id, cd, last)
?
Definition congen_distribution.f90:417
congen_distribution::wfn
subroutine, private wfn(nncsf, nadel, iidist, iidis3, nodi, ndi, cdi, ndel, pqnshl, x, ix, nx)
TODO...
Definition congen_distribution.f90:1816
congen_distribution::state
subroutine, private state(ns, x, last, nd, confpf)
?
Definition congen_distribution.f90:1604
congen_distribution::print5
subroutine, private print5(nd, ndi)
?
Definition congen_distribution.f90:1553
congen_distribution::getsa
subroutine, private getsa(ne, l, is, isz, m, nc, c, iso)
Form ne electrons in a shell coupled to (l,is).
Definition congen_distribution.f90:290
congen_distribution::packdet
subroutine, private packdet(iqn, ni, jqn, nj, ci, cj, nt)
Copy determinant data.
Definition congen_distribution.f90:1282
congen_distribution::print2
subroutine, private print2(iidis1)
?
Definition congen_distribution.f90:1410
congen_distribution::print4
subroutine, private print4(nd, x, ix)
?
Definition congen_distribution.f90:1509
congen_distribution::assign
subroutine, private assign(nshl, ndist, nused, refcon, excon, qnstor, nx, nftw)
Assign quantum numbers to real shells.
Definition congen_distribution.f90:54
congen_distribution::cgcoef
subroutine, private cgcoef(j1, j2, j3, m3, n, ms, c, intpfg)
Clebsch-Gordan coefficients.
Definition congen_distribution.f90:198
congen_pagecontrol
Page control.
Definition congen_pagecontrol.f90:26
congen_pagecontrol::taddl1
subroutine, public taddl1(lines, lt)
Check available lines.
Definition congen_pagecontrol.f90:193
congen_pagecontrol::newpg
subroutine, public newpg
Start new page of output.
Definition congen_pagecontrol.f90:121
congen_pagecontrol::space
subroutine, public space(lines)
Add blank lines to output.
Definition congen_pagecontrol.f90:144
congen_pagecontrol::addl
subroutine, public addl(lines)
Add blank lines to output.
Definition congen_pagecontrol.f90:50
congen_pagecontrol::taddl
subroutine, public taddl(lines, lt)
Add blank lines to output.
Definition congen_pagecontrol.f90:170