-
Notifications
You must be signed in to change notification settings - Fork 10
/
diatom.f90
17228 lines (17092 loc) · 643 KB
/
diatom.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module diatom_module
!
use accuracy
use timer
use functions, only : fanalytic_fieldT,define_complex_analytic_field
use symmetry, only : sym,SymmetryInitialize
use Lobatto, only : LobattoAbsWeights,derLobattoMat
use me_numer, only : ME_numerov
!
implicit none
! by Lorenzo Lodi
! Specifies variables needed for storing a temporary (scratch) copy
! of the input (useful for jumping around).
!
integer :: ierr
logical :: zDebug =.true. ! switching off this parameter will suppress many internal checks and should speed up calculations but make them more risky
character(len=wl) :: line_buffer
integer,parameter :: max_input_lines=500000 ! maximum length (in lines) of input. 500,000 lines is plenty..
! ! to avoid feeding in GB of data by mistake.
!
! main method used for solving the J=0 problem, i.e. to solve the one-dimentional Schrodinger equation.
! character(len=wl) :: solution_method = "5POINTDIFFERENCES" ! kinetic energy approximated by 5-point finite-diff. formula
character(len=wl) :: solution_method = "SINC" ! kinetic energy approximated by sinc DVR
!
! Type to describe different terms from the hamiltonian, e.g. potential energy, spin-orbit, <L^2>, <Lx>, <Ly> functions.
!
integer(ik),parameter :: verbose=5
integer(ik),parameter :: Nobjects = 33 ! number of different terms of the Hamiltonian
! (poten,spinorbit,L2,LxLy,spinspin,spinspino,bobrot,spinrot,diabatic,
! lambda-opq,lambda-p2q)
!
! In order to add a new field:
! 1. Change Nobjects
! 2. Add the name of the field to the CASE line where all fields are listed: case("SPIN-ORBIT","SPIN-ORBIT-X"....
! 3. Add the CASE-section describing the input of the new field
! 4. Add a similar CASE-section describing the input of the new field to the CASE("ABINITIO") part
! 5. Introduce a new name for the new field (object).
! 6. Define new array in type(fieldT),pointer .....
! 7. Add a new case to all cases where the fieds appear: at the end of the input subroutine,
! two times in in map_fields_onto_grid, in duo_j0, in sf_fitting
! 8. Add corresposponding check_and_print_coupling in map_fields_onto_grid
! 9. Add the corresponding name of the field to "use diatom_module,only" in "refinement"
!
! Current list of fields:
!
! case (1) poten(iterm)
! case (2) spinorbit(iterm)
! case (3) l2(iterm)
! case (4) lxly(iterm)
! case (5) spinspin(iterm)
! case (6) spinspino(iterm)
! case (7) bobrot(iterm)
! case (8) spinrot(iterm)
! case (9) diabatic(iterm)
! case (10) lambdaopq(iterm)
! case (11) lambdap2q(iterm)
! case (12) lambdaq(iterm)
! case (13) nac(iterm)
! case (14 - 20) reserved
! case (21) hfcc1(1) for Fermi contact, bF
! case (22) hfcc1(2) for nuclear spin - orbit, a
! case (23) hfcc1(3) for nuclear dipole - spin dipole, c
! case (24) hfcc1(4) for nuclear dipole - spind dipole, d
! case (25) hfcc1(5) for nuclear spin - rotaton, cI
! case (26) hfcc1(6) for electric dipole, eQq0
! case (27) hfcc1(7) for electric diople, eQq2
! case (28) reserved
! case(Nobjects-4) magnetictm(iterm)
! case(Nobjects-3) quadrupoletm(iterm)
! case(Nobjects-2) abinitio(iterm)
! case(Nobjects-1) brot(iterm)
! case(Nobjects) dipoletm(iterm)
!
!
! Lorenzo Lodi: it is necessary to repeat the character length specification in the array contructor for Fortran2003 conformance
character(len=wl),parameter :: CLASSNAMES(1:Nobjects) = (/ character(len=wl):: &
"POTEN","SPINORBIT", &
"L2", "L+", &
"SPIN-SPIN", &
"SPIN-SPIN-O", &
"BOBROT", &
"SPIN-ROT", &
"DIABATIC", &
"LAMBDAOPQ", "LAMBDAP2Q","LAMBDAQ", &
"NAC", &
"", "", "", "", "", "", "", & ! reserved
"HFCC-BF-1", "HFCC-A-1", "HFCC-C-1", "HFCC-D-1", "HFCC-CI-1", "HFCC-EQQ0-1", "HFCC-EQQ2-1", &
"", & ! reserved
"MAGNETIC",&
"QUADRUPOLE", &
"ABINITIO", &
"BROT", &
"DIPOLE"/)
!
! Lorenzo Lodi
! What follows is a bunch of temporary variables needed for computation of derivatives of pecs and other curves
real(kind=rk) :: fmmm,fmm,fm,f0,fp,fpp,fppp ! value of PEC on grid stencil
real(kind=rk) :: der1, der2, der3, der4 ! derivatives
real(kind=rk), parameter :: thresh_rel_err_min_x = 1e-13_rk ! threshold (relative error) for convergence of minimum of PEC
integer(kind=ik), parameter :: max_iter_min_search=20 ! maximum number of iterations for finding minimum of PEC
real(kind=rk), parameter :: h=2e-3_rk ! `h' is the step size in ang. for computing numerical derivatives of pecs
real(kind=rk) :: x0, x1 ! generic tmp variables
! Lorenzo Lodi -- text variables containing the symbols of the two atoms. The following formats are supported:
! 1) Chemical symbol. E.g. C, O, Na, Sc, Fe ....
! 2) Chemical symbol - atomic number, e.g.: C-14, O-16, Na-23
character(len=15) :: symbol1="Undefined", symbol2="Undefined"
!
!
type symmetryT
!
integer(ik) :: gu = 0
integer(ik) :: pm = 0
!
end type symmetryT
!
! type describing the parameter geneology
!
type linkT
!
integer(ik) :: iobject
integer(ik) :: ifield
integer(ik) :: iparam
!
end type linkT
!
! type describing the parameter fitting range
!
type rangeT
!
real(rk) :: min = -safe_max
real(rk) :: max = safe_max
logical :: set = .false.
!
end type rangeT
!
type braketT
integer(ik) :: ilambda = 0.0_rk
integer(ik) :: jlambda = 0.0_rk
real(rk) :: sigmai = 0.0_rk
real(rk) :: sigmaj = 0.0_rk
real(rk) :: value = 0.0_rk
end type braketT
!
! files with the eigenvectors
!
type eigenfileT
!
character(len=cl) :: dscr ! file with fingeprints and descriptions of each levels + energy values
character(len=cl) :: primitives ! file with the primitive quantum numbres
character(len=cl) :: vectors = 'eigen' ! eigenvectors stored here
!
end type eigenfileT
!
type weightT
character(cl) :: wtype = 'GRID'
real(rk) :: alpha = 0.001_rk
real(rk) :: Vtop = 10000.0_rk
end type weightT
!
integer(ik), parameter :: bad_value=-10000 ! default `impossible' value for some quantum numbers
!
type fieldT
!
character(len=cl) :: name='(unnamed)' ! Identifying name of the function (default to avoid garbled outputs)
character(len=cl) :: type='NONE' ! Identifying type of the function
character(len=cl) :: class='NONE' ! Identifying class of the function (poten, spinorbit,dipole,abinitio etc)
character(len=cl),pointer :: sub_type(:) ! Identifying type of sub-functions (used for coupled representations)
!
! variable used for GRID curves to specify interpolation and extrapolation options
!
character(len=cl) :: interpolation_type='CUBICSPLINES'
!
!character(len=cl) :: interpolation_type='QUINTICSPLINES'
!
integer(ik) :: iref ! reference number of the term as given in input (bra in case of the coupling)
integer(ik) :: jref ! reference number of the coupling term as given in input (ket in case of the coupling)
integer(ik) :: istate ! the actual state number (bra in case of the coupling)
integer(ik) :: jstate ! the actual state number (ket in case of the coupling)
character(len=cl) :: iTAG ! reference State TAG of the term as given in input (bra in case of the coupling), used to identify a state in the input
character(len=cl) :: jTAG ! reference State TAG of the term as given in input (ket in case of the coupling)
integer(ik) :: Nterms ! Number of terms or grid points
integer(ik) :: Nsubfunctions=0 ! Number of sub-functions
integer(ik),pointer :: Nsub_terms(:)=>null() ! Number of terms of sub-functions, used for coupled representations
integer(ik) :: Lambda =bad_value ! identification of the electronic state Lambda
integer(ik) :: Lambdaj=bad_value ! identification of the electronic state Lambda (ket)
integer(ik) :: multi ! identification of the electronic spin (bra for the coupling)
integer(ik) :: jmulti ! identification of the ket-multiplicity
real(rk) :: sigmai=real(bad_value,rk) ! the bra-projection of the spin
real(rk) :: sigmaj=real(bad_value,rk) ! the ket-projection of the spin
real(rk) :: spini ! electronic spin (bra component for couplings)
real(rk) :: spinj ! electronic spin of the ket vector for couplings
real(rk) :: omegai ! projection of J
real(rk) :: omegaj ! projection of J
integer(ik) :: iomega ! additional QN
integer(ik) :: jomega ! additional QN
integer(ik) :: ilevel ! additional QN
integer(ik) :: jlevel ! additional QN
complex(rk) :: complex_f=(1._rk,0._rk) ! defines if the term is imaginary or real
real(rk) :: factor=1.0_rk ! defines if the term is imaginary or real
real(rk) :: fit_factor=1.0 ! defines if the term is imaginary or real
real(rk),pointer :: value(:)=>null() ! Expansion parameter or grid values from the input
type(symmetryT) :: parity ! parity of the electronic state as defined by the molecular inversion (g,u),
! or laboratory inversion (+,-)
real(rk),pointer :: gridvalue(:) ! Expansion parameter or a grid value on the grid used inside the program
real(rk),pointer :: weight(:) ! fit (1) or no fit (0)
real(rk),pointer :: grid(:) ! grid value
real(rk),pointer :: matelem(:,:) ! matrix elements
real(rk) :: refvalue = 0 ! reference value will be used as a shift to be applied to the ab initio function used
! for the fitting constraints
character(len=cl),pointer :: forcename(:) ! The parameter name
integer(ik) :: nbrakets=0 ! total number of different combinations of lambda and sigma in matrix elements (maximum 4)
type(braketT) :: braket(4) ! here all possible combinations of values <+/-lambda,+/-sigma|A|+/-lambdaj,+/-sigma>
! can be listed
integer(ik) :: imin ! grid point index at which a potential energy curve is minimum
real(rk) :: Vimin ! V(imin), minimum of potential energy curve on the grid in cm-1
real(rk) :: rimin ! r_imin, r of the minimum of potential energy curve on the grid
integer(ik) :: imax ! grid point index at which a potential energy curve is maximum
real(rk) :: Vimax ! V(imax), maximum of potential energy curve on the grid in cm-1
logical :: zHasMinimum =.true. ! whether a potential energy curve has a minimum
real(rk) :: re ! pecs only: minimum for given pec, in angstroms.
real(rk) :: V0 ! pecs only: V(re) , in cm-1
real(rk) :: V1 ! pecs only: V'(re) , in cm-1 / angstroms (should be zero).
real(rk) :: V2 ! pecs only: V''(re) , in cm-1 / angstroms**2.
real(rk) :: V3 ! pecs only: V'''(re) , in cm-1 / angstroms**3.
real(rk) :: V4 ! pecs only: V''''(re), in cm-1 / angstroms**4.
real(rk) :: we ! pecs only: harmonic vibrational freq, cm-1.
real(rk) :: B0 ! pecs only: v=0 vibrational constant vibrational, cm-1.
real(rk) :: xe ! pecs only: anharmonicity constant (pure number)
real(rk) :: alphae ! pecs only: Coriolis coupling constant, cm-1.
real(rk) :: Debar ! pecs only: Centrifugal distortion constant, cm-1.
real(rk) :: Y00 ! pecs only: ZPE correction, cm-1.
real(rk) :: Omega_min ! pecs only: minimum physically possible value for |Omega|=|Lambda+Sigma|
real(rk) :: approxEJ0 ! pecs only: approximate J=0, v=0 energy (no couplings)
real(rk) :: approxEJmin ! pecs only: approximate J=Jmin, v=0 energy (no couplings)
procedure (fanalytic_fieldT),pointer, nopass :: fanalytic_field => null()
!
type(linkT),pointer :: link(:) ! address to link with the fitting parameter in a different object in the fit
logical :: morphing = .false. ! When morphing is the field is used to morph the ab initio couterpart
! towards the final object
logical :: molpro = .false. ! The object is given in the molpro representaion (|x> and |y>)
integer(ik) :: ix_lz_y = 1000 ! The value of the matrix element (-I)<x|Lz|y> for i-state, 1000 is for undefined
integer(ik) :: jx_lz_y = 1000 ! The value of the matrix element (-I)<x|Lz|y> for j-state, 1000 is for undefined
type(weightT) :: weighting ! When morphing is the field is used to morph the ab initio couterpart
!
character(len=cl) :: integration_method='DVR-SINC' ! Identifying the type of the integration method used specifically for this field, DVR-SINC default
!
real(rk) :: adjust_val = 0.0_rk
real(rk) :: asymptote = 0.0_rk ! reference asymptote energy used e.g. in the renormalisation of the continuum wavefunctions to sin(kr)
logical :: adjust = .false. ! Add constant adjust_val to all fields
integer(ik) :: iabi = 0 ! abinitio field's id if associated
type(rangeT),pointer :: fit_range(:) ! fitting range, min..max
!
end type fieldT
!
type FieldListT
TYPE(fieldT), POINTER :: field(:)
INTEGER :: num_field = 0
PROCEDURE(), POINTER, NOPASS :: hfcc_matrix_element => NULL()
end type FieldListT
type jobT
logical :: select_gamma(4) ! the diagonalization will be done only for selected gamma's
integer(ik) :: nroots(4)=1e9_rk ! number of the roots to be found in variational diagonalization with syevr
integer(ik) :: maxiter = 1000 ! maximal number of iterations in arpack
real(rk) :: tolerance = 0.0_rk ! tolerance for arpack diagonalization, 0 means the machine accuracy
real(rk) :: upper_ener = 1e9_rk ! upper energy limit for the eigenvalues found by diagonalization with syevr
real(rk) :: thresh = -1e-18_rk ! thresh of general use
real(rk) :: zpe=0.0_rk ! zero-point-energy
logical :: shift_to_zpe = .true. !
character(len=cl) :: diagonalizer = 'SYEV'
character(len=cl) :: molecule = 'H2'
character(len=cl) :: contraction = 'VIB' ! contraction
real(rk),pointer :: vibenermax(:) ! contraction parameter: energy
integer(ik),pointer :: vibmax(:) ! contraction parameter: vibration quantum
real(rk) :: potmin ! absolute minimum of the PEC with the lowest rotational-vibrational state
logical :: zShiftPECsToZero = .true.
logical :: zEchoInput = .true.
logical :: zExclude_JS_coupling =.false. ! If set to true will disable J.S coupling (aka S-uncoupling)
integer(ik) :: total_parameters =0 ! total number of parameters used to define different hamiltonian fields
real(rk) :: degen_threshold = 1e-6_rk
real(rk),pointer :: j_list(:) ! J values processed
integer(ik) :: nJ = 1 ! Number of J values processed
character(len=cl) :: IO_eigen = 'NONE' ! we can either SAVE to or READ from the eigenfunctions from an external file
character(len=cl) :: IO_density = 'NONE' ! we can SAVE reduced density to an external file
character(len=cl) :: IO_dipole = 'NONE' ! we can either SAVE to or READ from an external file the dipole moment
character(len=cl) :: basis_set = 'NONE' ! we can keep the vibrational basis functions for further usage
! matrix elements on the contr. basis
type(eigenfileT) :: eigenfile
character(len=cl) :: symmetry = 'CS(M)' ! molecular symmetry
real(rk) :: diag_L2_fact = 1._rk ! specifies the convention used for the diagonal contribution
! due to L^2 = LxLx+LyLy+LzLz
logical,pointer :: isym_do(:) ! process or not the symmetry in question
logical :: intensity ! switch on the intensity calculations
logical :: print_vibrational_energies_to_file = .false. ! if .true. prints to file
logical :: print_rovibronic_energies_to_file = .false. ! if .true. prints to file
logical :: print_pecs_and_couplings_to_file = .false. ! if .true. prints to file
logical :: assign_v_by_count = .true.
logical :: legacy_version = .false.
!
end type jobT
!
type gridT
integer(ik) :: npoints = 1000 ! grid size
real(rk) :: rmin = 1.0,rmax=3.00 ! range of the grid
real(rk) :: step = 1e-2 ! step size
real(rk) :: alpha = 0.2 ! grid parameter
real(rk) :: re = 1.0 ! grid parameter
integer(ik) :: nsub = 0 ! grid type parameter (0=uniformly spaced)
real(rk),pointer :: r(:)=>null() ! the molecular geometry at the grid point
end type gridT
type quantaT
real(rk) :: I1 ! nuclear spin 1
real(rk) :: F1 ! \hat{F1} = \hat{I1} + \hat{J}
INTEGER(ik) :: index_F1 ! index of the F1 in F1_list
real(rk) :: Jrot ! J - real
integer(ik) :: irot ! index of the J value in J_list
integer(ik) :: istate ! e-state
integer(ik) :: imulti ! imulti = 1,..,multiplicity
real(rk) :: sigma ! spin-projection = -spin,-spin+1,..,spin-1,spin
real(rk) :: omega ! sigma+lambda
real(rk) :: spin ! spin
integer(ik) :: ilambda ! ilambda
integer(ik) :: v = 0 ! vibrational quantum
integer(ik) :: ivib = 1 ! vibrational quantum number counting all vib-contracted states
integer(ik) :: ilevel = 1 ! primitive quantum
integer(ik) :: iroot ! running number
integer(ik) :: iF1_ID ! running number within the same F and parity
integer(ik) :: iJ_ID ! running number within the same J
integer(ik) :: iparity = 0
integer(ik) :: igamma = 1
integer(ik) :: iomega = 1 ! countig number of omega
character(len=cl) :: name ! Identifying name of the function
logical :: bound = .true. ! is this state bound or unbound
character(len=cl) :: iTAG ! reference State TAG of the term as given in input (bra in case of the coupling), used to identify a state in the input
end type quantaT
!
type eigenT
real(rk),pointer :: vect(:,:) ! the eigenvector in the J=0 contracted representaion
real(rk),pointer :: val(:) ! the eigenvalue
type(quantaT),pointer ::quanta(:) ! the quantum numbers
integer(ik) :: Nlevels
integer(ik) :: Ndimen
integer(ik) :: Nbound
end type eigenT
!
type basisT
type(quantaT),pointer :: icontr(:) ! the quantum numbers
integer(ik) :: Ndimen
end type basisT
!
type actionT
!
logical :: fitting = .false.
logical :: intensity = .false.
logical :: frequency = .false.
logical :: matelem = .false.
logical :: raman = .false.
logical :: magdipole = .false.
logical :: quadrupole = .false.
logical :: RWF = .false.
logical :: overlap = .false.
logical :: hyperfine = .false.
logical :: save_eigen_J = .false.
!
end type actionT
!
type obsT
!
real(rk) :: Jrot
real(rk) :: Jrot_ ! J - real (lower)
integer(ik) :: irot ! index of the J value in J_list
integer(ik) :: irot_
integer(ik) :: symmetry
integer(ik) :: N
integer(ik) :: N_ ! N lower
integer(ik) :: iparity
integer(ik) :: iparity_ ! parity +/- (0,1) of the lower state lower
real(rk) :: energy
real(rk) :: frequency
real(rk) :: weight
type(quantaT) :: quanta
type(quantaT) :: quanta_ ! quantum numbers of the lower state
!
end type obsT
type thresholdsT
real(rk) :: intensity = -1e0 ! threshold defining the output intensities
real(rk) :: linestrength = -1e0 ! threshold defining the output linestrength
real(rk) :: dipole = -1e0 ! threshold defining the output linestrength
real(rk) :: quadrupole = -1e0 ! threshold defining the output linestrength
real(rk) :: coeff = -1e0 ! threshold defining the eigenfunction coefficients
! taken into account in the matrix elements evaluation.
real(rk) :: bound_density = sqrt(small_) ! threshold defining the unbound state density
! calculated at the edge of the box whioch must be small for bound states
real(rk) :: bound_aver_density = sqrt(small_) ! is bound_density/deltaR_dens
real(rk) :: deltaR_dens = 0.5_rk ! small interval for computing the state density (Angstrom)
!
end type thresholdsT
!
type landeT
real, allocatable :: ar(:)
end type landeT
!
type IntensityT
logical :: do = .false. ! process (.true.) or not (.false.) the intensity (or TM) calculations
character(cl) :: action ! type of the intensity calculations:
! absorption, emission, tm (transition moments),
! raman, and so on.
real(rk) :: temperature ! temperature in K
real(rk) :: part_func=0 ! partition function
real(rk) :: ZPE=0 ! zero point energy
type(thresholdsT) :: threshold ! different thresholds
real(rk),pointer :: gns(:) ! nuclear stat. weights
integer(ik),pointer :: isym_pairs(:) ! numbers defining symmetry pairs with allowed transitions, analogous to gns
real(rk) :: freq_window(1:2) ! frequency window (1/cm)
real(rk) :: erange_low(1:2) ! energy range for the lower state
real(rk) :: erange_upp(1:2) ! energy range for the upper state
real(rk) :: J(1:2) ! range of J-values, from..to; in order to avoid double counting of transitions
! in the calculations it is always assumed that
! J1<=J_lower<=J2 and J1<=J_upper<J2;
!
type(quantaT) :: lower ! lower state range of the quantun numbers employed
type(quantaT) :: upper ! upper state range of the quantun numbers employed
! in intensity calculations; (imode,1:2),
! where 1 stands for the beginning and 2 for the end.
!
integer(ik) :: swap_size = 0 ! the number of vectors to keep in memory
character(cl) :: swap = "NONE" ! whether save the compacted vectors or read
character(cl) :: swap_file ="compress" ! where swap the compacted eigenvectors to
character(cl) :: linelist_file="NONE" ! filename for the line list (filename.states and filename.trans)
integer(ik) :: int_increm = 1e9 ! used to print out the lower energies needed to select int_increm intensities
real(rk) :: factor = 1.0d0 ! factor <1 to be applied the maxsize of the vector adn thus to be shrunk
logical :: matelem =.false. ! switch for the line-strenth-type matelems (matrix elements of the dipole moment)
logical :: lande_calc = .false. ! checks whether calculation for Lande should be conducted
logical :: overlap = .false. ! print out overlap integrals (Franck-Condon)
logical :: tdm = .true. ! print out dipole transition moments
logical :: tqm = .true. ! print out quadrupole transition moments
integer(ik) :: Npoints = 10 ! used for cross sections grids
real(rk) :: gamma = 0.05_rk ! Lorentzian FWHM, needed for cross-sections
integer(ik) :: N_RWF_order = 1 ! Expansion order of the matrix fraction needed for RWF
character(cl) :: RWF_type="GAUSSIAN" ! Type of RWH
logical :: renorm = .false. ! renormalize the continuum/unbound wavefunctions to sin(kr) for r -> infty
logical :: bound = .false. ! filter bound states
logical :: unbound = .false. ! filter and process unbound upper states only
logical :: interpolate = .false. ! use interpolation of Einstein coefficients on a grid of frequencies
logical :: states_only = .false. ! Only .states file is generated while .trans is skipped. Equivalent to setting negative freq-window
real(rk) :: asymptote=safe_max ! Lowest asymptote
!
end type IntensityT
!
type matrixT
!
real(rk),pointer :: matrix(:,:)
integer(ik),pointer :: irec(:)
!
end type matrixT
!
type fittingT
!
logical :: run
integer(ik) :: nJ = 1 ! Number of J values processed
real(rk),pointer :: j_list(:) ! J-values processed in the fit
integer(ik) :: iparam(1:2) = (/1,100000/)
integer(ik) :: itermax = 500
integer(ik) :: Nenergies = 1
integer(ik) :: parmax =0 ! total number of all parameters used
real(rk) :: factor = 1.0_rk
real(rk) :: target_rms = 1e-8
real(rk) :: robust = 0
character(len=cl) :: geom_file = 'pot.fit'
character(len=cl) :: output_file = 'fitting'
character(len=cl) :: fit_type = 'LINUR' ! to switch between fitting methods.
real(rk) :: threshold_coeff = -1e-18
real(rk) :: threshold_lock = -1e-18
real(rk) :: threshold_obs_calc = -1e-16
real(rk) :: zpe=0
logical :: shift_to_zpe = .true. !
real(rk) :: fit_scaling=1.0_rk ! scaling the fitting correction with this factor >0 and <1
integer(ik) :: linear_search = 0 ! use linear scaling to achieve better convergence with the Armijo condition
type(obsT),pointer :: obs(:) ! experimental data
!
!type(paramT),pointer :: param(:) ! fitting parameters
!
end type fittingT
!
type fieldmapT
integer(ik) :: Nfields
end type fieldmapT
!
! This type is designed for the Omega-state representations
type Omega_gridT
integer(ik) :: Nstates
real(rk) :: omega
real(rk),pointer :: energy(:,:)=>null() ! the Omega-potential energy at each grid point
real(rk),pointer :: vector(:,:,:)=>null() ! the Omega-vector at each grid point in the lambda-sigma space
type(quantaT),pointer :: QN(:)=>null() ! quantum numbers
type(quantaT),pointer :: basis(:)=>null() ! basis set
end type Omega_gridT
!
! This type is for the vibronic contacted eigensolutions in Omega represetation
type contract_solT
integer(ik) :: Ndimen
real(rk),pointer :: Energy(:)
real(rk),pointer :: vector(:,:)=>null() ! the Omega-vector at each grid point in the lambda-sigma space
integer(ik),pointer :: ilevel(:)=>null() ! level index
end type contract_solT
type F1_hyperfine_steup_T
real(rk) :: I1
real(rk) :: fit_bound_ratio
character(cl) :: fit_optimization_algorithm
end type F1_hyperfine_steup_T
!
integer, parameter :: trk = selected_real_kind(12)
integer,parameter :: jlist_max = 500
integer,parameter :: states_max = 50
type(fieldT),pointer :: poten(:)=>null(),spinorbit(:)=>null(),l2(:)=>null(),lxly(:)=>null(),abinitio(:)=>null(),&
dipoletm(:)=>null(),spinspin(:)=>null(),spinspino(:)=>null(),bobrot(:)=>null(),spinrot(:)=>null(),&
diabatic(:)=>null(),lambdaopq(:)=>null(),lambdap2q(:)=>null(),lambdaq(:)=>null(),nac(:)=>null()
type(fieldT),pointer :: brot(:)=>null(),quadrupoletm(:)=>null(),magnetictm(:)=>null()
!
! Fields in the Omega representation
!
integer(ik) :: Nspins,NLplus_omega,NSplus_omega,NSR_omega,NBob_omega,Nomegas,Np2q_omega,Nq_omega,NKin_omega,NBRot_omega,NNAC_omega
integer(ik) :: NDiab_omega,NDipole_omega
type(Omega_gridT),allocatable :: omega_grid(:)
integer(ik),allocatable :: iLplus_omega(:,:,:,:),iSplus_omega(:,:,:,:),iSR_omega(:,:,:,:),iBOB_omega(:,:,:)
integer(ik),allocatable :: iP2Q_omega(:,:,:,:),iQ_omega(:,:,:,:),iKin_omega(:,:,:),iBRot_omega(:,:,:),iNAC_omega(:,:,:),&
iDiab_omega(:,:,:),iDipole_omega(:,:,:,:)
!
type(fieldT),pointer :: l_omega_obj(:)=>null(),s_omega_obj(:)=>null(),sr_omega_obj(:)=>null(),bob_omega_obj(:)=>null()
type(fieldT),pointer :: p2q_omega_obj(:)=>null(),q_omega_obj(:)=>null(),kin_omega_obj(:)=>null(),brot_omega_obj(:)=>null()
type(fieldT),pointer :: nac_omega_obj(:)=>null(),Diab_omega_obj(:)=>null(),Dipole_omega_obj(:)=>null()
!
type(jobT) :: job
type(gridT) :: grid
type(quantaT),allocatable :: quanta(:)
integer(ik),allocatable :: iquanta2ilevel(:,:,:)
real(rk),allocatable :: r(:), d2dr(:), r2sc(:),z(:)
type(actionT) :: action ! defines different actions to perform
type(fittingT) :: fitting
type(IntensityT) :: Intensity
!type(symmetryT) :: sym
!
integer(ik) :: nestates=-1,Nspinorbits,Ndipoles,Nlxly,Nl2,Nabi,Ntotalfields=0,Nss,Nsso,Nbobrot,Nsr,Ndiabatic,&
Nlambdaopq,Nlambdap2q,Nlambdaq,Nnac,vmax,nQuadrupoles,NBrot,nrefstates = 1,nMagneticDipoles
real(rk) :: m1=-1._rk,m2=-1._rk ! impossible, negative initial values for the atom masses
real(rk) :: jmin,jmax,amass,hstep,Nspin1=-1.0,Nspin2=-1.0
real(rk) :: jmin_global
!
!type(fieldT),pointer :: refined(:)
type(fieldmapT) :: fieldmap(Nobjects)
type(eigenT),allocatable :: eigen(:,:)
!
type(basisT),allocatable :: basis(:)
!
logical :: gridvalue_allocated = .false.
logical :: fields_allocated = .false.
real(rk),parameter :: enermax = safe_max ! largest energy allowed
real(rk), allocatable :: vibrational_contrfunc(:,:)
type(quantaT), allocatable :: vibrational_quantum_number(:)
integer(ik) :: vibrational_totalroots
type(F1_hyperfine_steup_T) :: F1_hyperfine_setup
integer(ik), parameter :: GLOBAL_NUM_HFCC_OBJECT = 7
type(FieldListT) :: hfcc1(GLOBAL_NUM_HFCC_OBJECT)
!
real(rk),allocatable :: overlap_matelem(:,:)
!
public ReadInput,poten,spinorbit,l2,lxly,abinitio,brot,map_fields_onto_grid,fitting,&
jmin,jmax,vmax,fieldmap,Intensity,eigen,basis,Ndipoles,dipoletm,linkT,rangeT,three_j,quadrupoletm,&
magnetictm,l_omega_obj,s_omega_obj,sr_omega_obj,brot_omega_obj,p2q_omega_obj,q_omega_obj,kin_omega_obj,&
nac_omega_obj,overlap_matelem,Diab_omega_obj,Dipole_omega_obj
!
save grid, Intensity, fitting, action, job, gridvalue_allocated, fields_allocated, hfcc1
!
contains
!
subroutine ReadInput
!
use input
!
integer(ik) :: iobject(Nobjects)
integer(ik) :: ipot=0,iso=0,ncouples=0,il2=0,ilxly=0,iabi=0,idip=0,iss=0,isso=0,isr=0,idiab=0,iquad=0,imagnetic=0
integer(ik) :: Nparam,alloc,iparam,i,j,iobs,i_t,iref,jref,istate,jstate,istate_,jstate_,item_,ibraket,iabi_,&
iobj,iclass_,ielement,nstate_listed
integer(ik) :: Nparam_check !number of parameters as determined automatically by duo (Nparam is specified in input).
logical :: zNparam_defined ! true if Nparam is in the input, false otherwise..
integer(ik) :: itau,lambda_,x_lz_y_,iobject_,inac
logical :: integer_spin = .false., matchfound
real(rk) :: unit_field = 1.0_rk,unit_adjust = 1.0_rk, unit_r = 1.0_rk,spin_,jrot2
real(rk) :: f_t,jrot,j_list_(1:jlist_max)=-1.0_rk,omega_,sigma_,hstep = -1.0_rk
real(rk) :: bound_density
!
character(len=cl) :: w,ioname,iTAG,jTAG,States_list(states_max)
character(len=wl) :: large_fmt
!
integer(ik) :: iut ! iut is a unit number.
!
type(fieldT),pointer :: field,field_
logical :: eof,include_state,allgrids
logical :: symmetry_defined=.false.,skip_diabatic,skip_magnetic
integer :: ic,ierr
!
! -----------------------------------------------------------
!
!
! read the general input
! by Lorenzo Lodi
! read everything from std input and write to a temporary (scratch) file.
!
write(ioname, '(a, i4)') 'write to a temporary (scratch) file.'
call IOstart(trim(ioname), iut)
!
open(unit=iut, status='scratch', action='readwrite')
write(large_fmt, '(A,i0,A)') '(A', wl, ')'
trans_loop: do i=1, max_input_lines
read(unit=*,fmt=large_fmt,iostat=ierr) line_buffer
if(ierr /=0) exit trans_loop
write(iut, '(a)') trim(line_buffer)
! This is a hack; I need to know if to echo the input or not before processing it
! The option 'do_not_echo_input' is dealt with here as a special case
line_buffer = adjustl(line_buffer) ! remove leading spaces
do j=1, len(trim(line_buffer)) ! convert to uppercase
ic = ichar( line_buffer(j:j))
if( ic >= 97) line_buffer(j:j) = achar(ic-32)
enddo
if( trim(line_buffer) == 'DO_NOT_ECHO_INPUT' ) job%zEchoInput = .false.
enddo trans_loop
rewind(iut)
!
!
! default constants
!
jmin = 0 ; jmax = 0
!
! To count objects
iobject = 0
!
nstate_listed = 0
!
if( job%zEchoInput) then
write(out,"('(Transcript of the input --->)')")
call input_options(echo_lines=.true.,error_flag=1)
else
call input_options(echo_lines=.false.,error_flag=1)
endif
do
zNparam_defined = .false. ! For each input block, set number of points/params to undefined
call read_line(eof,iut) ; if (eof) exit
call readu(w)
select case(w)
!
case("STOP","FINISH","END")
!
exit
!
case("") ! do nothing in case of blank lines
!
case("PRINT_PECS_AND_COUPLINGS_TO_FILE")
job%print_pecs_and_couplings_to_file = .true.
!
case("PRINT_VIBRATIONAL_ENERGIES_TO_FILE")
job%print_vibrational_energies_to_file = .true.
!
case("PRINT_ROVIBRONIC_ENERGIES_TO_FILE")
!
job%print_rovibronic_energies_to_file = .true.
case("DO_NOT_ECHO_INPUT")
!
job%zEchoInput = .false.
case("DO_NOT_SHIFT_PECS")
!
job%zShiftPECsToZero = .false.
case("DO_NOT_INCLUDE_JS_COUPLING")
!
job%zExclude_JS_coupling = .true.
!
case("ASSIGN_V_BY_COUNT")
!
job%assign_v_by_count = .true.
!
case("ASSIGN_V_BY_CONTRIBUTION")
!
job%assign_v_by_count = .false.
!
case("LEGACY","OLD-VERSION","VERSION")
!
job%legacy_version = .true.
!
if (item>1) then
call readi(item_)
if (item_>2018) job%legacy_version = .false.
endif
!
case ("SOLUTIONMETHOD")
!
call readu(w)
solution_method = trim(w)
!
case ("L2CONVENTION")
!
call readu(w)
!
select case(w)
case ("SPECIFY_L^2","SPECIFY_L**2" ,"DEFAULT")
job%diag_L2_fact = 1._rk !default case
case ("SPECIFY_LX^2_PLUS_LY^2", "SPECIFY_LX**2_PLUS_LY**2" )
job%diag_L2_fact = 0._rk !Lorenzo's choice
case default
call report ("Unrecognized L2CONVENTION "//trim(w)// ". " // &
"Implemented: DEFAULT, SPECIFY_L**2, SPECIFY_L^2, SPECIFY_LX^2_PLUS_LY^2" // &
", SPECIFY_LX**2_PLUS_LY**2")
end select
!
case ("MEM","MEMORY")
!
call readf(memory_limit)
!
if (nitems<2) call report("Please specify the memory units",.true.)
!
call readu(w)
!
select case(w)
!
case default
!
call report("Unexpected argument in MEMORY",.true.)
!
case("TB","T")
!
memory_limit = memory_limit*1000.0_rk**2
!
case("GB","G")
!
memory_limit = memory_limit*1000.0_rk
!
case("MB","M")
!
memory_limit = memory_limit
!
case("KB","K")
!
memory_limit = memory_limit/1000.0_rk
!
case("B")
!
memory_limit = memory_limit/1000.0_rk**2
!
end select
!
case ("ATOMS") ! chemical symbols of the atoms (possibly with atomic numbers)
!
call reada(symbol1)
call reada(symbol2)
!
case ("MASSES","MASS")
!
call readf(m1)
call readf(m2)
!
case ("MOLECULE","MOL")
!
call readu(w)
!
select case(w)
!
case ("C2","ALO","X2","XY","CO","CAO","NIH","MGH")
!
job%molecule = trim(w)
!
case default
!
! write (out,"(' I see this molecule for the first time.')")
!
!call report ("Unrecognized unit name "//trim(w)//"implemented: (C2,ALO,X2, XY)",.true.)
!
end select
!
case('J_LIST','JLIST','JROT','J')
!
jmin = 1e6
jmax = -1.0
integer_spin = .false.
i = 0
do while (item<Nitems.and.i<jlist_max)
!
i = i + 1
!
call readu(w)
!
if (trim(w)/="-") then
!
read(w,*) jrot
!
j_list_(i) = jrot
!
else
!
call readf(jrot2)
!
do while (item<=Nitems.and.nint(2.0*jrot)<nint(2.0*jrot2))
!
jrot = jrot + 1.0
j_list_(i) = jrot
i = i + 1
!
enddo
i = i - 1
!
endif
!
if (i==1.and.mod(nint(2.0_rk*jrot+1.0_rk),2)==1) integer_spin = .true.
!
if (i>1.and.mod(nint(2.0_rk*jrot+1.0_rk),2)==0.and.integer_spin) then
!
call report("The multiplicities of J-s in J-list/Jrot are inconsistent",.true.)
!
endif
!
jmin = min(jmin,j_list_(i))
jmax = max(jmax,j_list_(i))
!
enddo
!
job%nJ = i
!
allocate(job%j_list(i),stat=alloc)
!
job%J_list(1:i) = J_list_(1:i)
!
!case ("JROT")
! !
! call readf(jmin)
! !
! if (nitems>2) then
! !
! call readf(jmax)
! !
! else
! !
! jmax = jmin
! !
! endif
! !
! ! check the multiplicity
! !
! if (mod(nint(2.0_rk*jmin+1.0_rk),2)==1) integer_spin = .true.
! !
! if (mod(nint(2.0_rk*jmax+1.0_rk),2)==0.and.integer_spin) then
! !
! call report("The multiplicities of jmin and jmax are inconsistent",.true.)
! !
! endif
! !
case ("NSTATES","NESTATES","STATES")
!
if (nestates>-1) then
!
call report("Conflicting input: Nstates and States cannot be used at the same time: ",.true.)
!
endif
!
select case (w)
!
case ("STATES")
!
i = 0
!
do while (item<Nitems.and.i<states_max)
!
i = i + 1
call reada(States_list(i))
!
enddo
!
nestates = i
!
nstate_listed = i
!
case default
!
call readi(nestates)
!
end select
!
if (nestates<1) call report("nestates cannot be 0 or negative",.true.)
!
! the maximum number of couplings possible, assuming each state may be doubly degenerate
!
ncouples = 2*nestates*(2*nestates-1)/2
!
! allocate all fields representing the hamiltonian: PEC, spin-orbit, L2, LxLy etc.
!
allocate(poten(nestates),spinorbit(ncouples),l2(ncouples),lxly(ncouples),spinspin(ncouples),spinspino(ncouples), &
bobrot(nestates),spinrot(nestates),job%vibmax(nestates),job%vibenermax(nestates),diabatic(ncouples),&
lambdaopq(nestates),lambdap2q(nestates),lambdaq(nestates),nac(ncouples),quadrupoletm(ncouples),&
magnetictm(ncouples),stat=alloc)
do i = 1, GLOBAL_NUM_HFCC_OBJECT
allocate(hfcc1(i)%field(nestates), stat=alloc)
end do
!
! initializing the fields
!
job%vibmax = 1e8
job%vibenermax = enermax
!
allocate(abinitio(nestates*Nobjects+5*ncouples),stat=alloc)
!
case ("NREFSTATES")
!
call readi(nrefstates)
!
if (nrefstates>nestates) call report("nrefstates cannot be larger than nestates",.true.)
!
case ("GRID")
!
call read_line(eof,iut) ; if (eof) exit
call readu(w)
!
do while (trim(w)/="".and.trim(w)/="END")
!
select case(w)
!
case ("NPOINTS")
!
call readi(grid%npoints)
!
case ("STEP")
!
call readf(hstep)
!
case ("NSUB","TYPE")
!
call readi(grid%nsub)
!
case ("RE","REF")
!
call readf(grid%re)
!
case ("ALPHA")
!
call readf(grid%alpha)
!
case("RANGE")
!
call readf(grid%rmin)
call readf(grid%rmax)
!
case default
!
call report ("Unrecognized keyword in GRID: "//trim(w),.true.)
!
end select
!
if (hstep>0.and.grid%npoints/=0) then
write(out,"('Illegal grid-input: npoints and step should not appear together')")
stop "Illegal grid-input: npoints and step should not appear together"
endif
!
if (hstep>0) then
grid%npoints = (grid%rmax - grid%rmin)/hstep+1
endif
!
call read_line(eof,iut) ; if (eof) exit
call readu(w)
!
enddo
!
if (trim(w)/="".and.trim(w)/="END") then
!
call report ("Unrecognized keyword in GRID: "//trim(w),.true.)
!
endif
!
case ("CONTRACTION","VIBRATIONS","VIBRATIONALBASIS")
!
call read_line(eof,iut) ; if (eof) exit
call readu(w)
!
do while (trim(w)/="".and.trim(w)/="END")
!
select case(w)
!
case('VIB','ROT','OMEGA')
!
job%contraction = trim(w)
!
case ("VMAX","VIBMAX","NMAX")
!
vmax = 1
!
if (Nitems<2) job%vibmax = grid%npoints-1