-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathsurfacehop.F90
1316 lines (1079 loc) · 48.4 KB
/
surfacehop.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
! Driver routines for Surface Hopping dynamics
! D. Hollas, O. Svoboda, P. Slavicek, M. Oncak
module mod_sh
use mod_const, only: DP
use mod_files, only: stdout, stderr
use mod_error, only: fatal_error
use mod_sh_integ
implicit none
private
! INPUT PARAMETERS
public :: read_sh_input, print_sh_input
public :: inac
! Main subroutines
public :: surfacehop, sh_init
! Helper subroutines
public :: move_vars, get_nacm, write_nacmrest, read_nacmrest
public :: check_CIVector
public :: set_current_state
! Variables holding SH state
! TODO: Make these protected and write Set methods
public :: en_array, tocalc
public :: nacx, nacy, nacz
! Initial electronic state
integer :: istate_init = -1
! Numerical integrator for SH wavefunction
! (default is Butcher 5-th order)
character(len=50) :: integ = 'butcher'
! Number of substeps for integrating electronic Schrodinger eq.
integer :: substep = 100
! Controls calculations of Non-adiabatic Couplings (NAC)
! couplings = 'analytic' (inac=0) - Analytical NAC (default)
! couplings = 'baeck-an' (inac=1) - Baeck-An couplings
! couplings = 'none' (inac=2) - Do not compute couplings
integer :: inac = 0 ! for working within the code
character(len=50) :: couplings = 'analytic' ! for reading the input file
! energy history array and time-derivate couplings (sigma_ba) necessary for Beack-An couplings
real(DP), allocatable :: en_hist_array(:, :), sigma_ba(:, :) ! sigma_ba is actually the dotproduct
! 1 - Turn OFF hopping
integer :: nohop = 0
! How to adjust velocity after hop:
! velocity_rescaling = 'nac_then_velocity' (adjmom=0) - Adjust velocity along the NAC vector, if not possible,
! try the velocity vector (default)
! velocity_rescaling = 'velocity' (adjmom=1) - Rescale along the velocity vector
integer :: adjmom = 0 ! for working within the code
character(len=50) :: velocity_rescaling = 'nac_then_velocity' ! for reading the input file
! 1 - Reverse momentum direction after frustrated hop
integer :: revmom = 0
! Numerical accuracy of MOLPRO NACME, 10^-nac_accu
! Default accuracy (equals MOLPRO defaults)
integer :: nac_accu1 = 7
! If calculation fails with default, we retry with nac_accu2
integer :: nac_accu2 = 5
! Decoherence correction parameter (a.u.)
! To turn off decoherence correction, set decoh_alpha = 0.0D0
real(DP) :: decoh_alpha = 0.1D0
! Compute NACME only if states are less then deltaE apart (eV)
! The default value (5 eV) is very conservative.
real(DP) :: deltaE = 5.0D0
! Compute NACME only if at least one of the two states has population > popthr
real(DP) :: popthr = 0.001D0
! Thresholds for energy conservations (in eV)
! The default values are too permisive
! you should definitely tighten them in production simulations!
! Stop simulation if total energy difference in successive steps is > energydifthr
real(DP) :: energydifthr = 1.0D0
! Stop simulation if total energy drift since the beginning exceeds energydriftthr
real(DP) :: energydriftthr = 1.0D0
! Surface Hopping with Induced Transitions,
! typically used with methods that can't describe S1S0 intersections,
! such as TDDFT or ADC2.
! In this case, we force a hop to S0 when the energy difference
! between S1 and S0 drops below user defined threshold.
logical :: SHwIT = .false.
real(DP) :: dE_S0S1_thr = 0.0D0 !eV
! NA Couplings
real(DP), allocatable :: nacx(:, :, :), nacy(:, :, :), nacz(:, :, :)
! *old variables holds data from the previous step
real(DP), allocatable :: nacx_old(:, :, :), nacy_old(:, :, :), nacz_old(:, :, :)
real(DP), allocatable :: en_array(:), en_array_old(:)
! Initial absolute electronic energy, needed for monitoring energy drift
real(DP) :: entot0
! nstate x nstate matrix to determine which derivatives to compute
! off-diagonal elements correspond to NAC
! diagonal elements correspond to electronic gradients
integer, allocatable :: tocalc(:, :)
! Current electronic state
integer, public, protected :: istate
! Do not consider hopping into this state
! Useful e.g. for ethylene 2-state SA3 dynamics
! This is a bit of an ugly hack, it would be more general to have an array
! of states that are calculated but ignored.
integer :: ignore_state = 0
namelist /sh/ istate_init, nstate, substep, deltae, integ, couplings, nohop, phase, decoh_alpha, popthr, ignore_state, &
nac_accu1, nac_accu2, popsumthr, energydifthr, energydriftthr, velocity_rescaling, revmom, &
shwit, dE_S0S1_thr, correct_decoherence
save
contains
! Initialization routine (ugly, but works)
subroutine sh_init(x, y, z, vx, vy, vz)
use mod_const, only: AUtoEV
use mod_general, only: irest, natom, pot
use mod_interfaces, only: force_clas
use mod_kinetic, only: ekin_v
use mod_files, only: nacmefile_init
real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :)
real(DP) :: dum_fx(size(x, 1), size(x, 2))
real(DP) :: dum_fy(size(y, 1), size(y, 2))
real(DP) :: dum_fz(size(z, 1), size(z, 2))
real(DP) :: dum_eclas
write (stdout, *) ''
write (stdout, '(A)') 'Initializing Surface Hopping'
call check_sh_parameters()
! TODO: Write a print_sh_header()
! with all major parameters included and explained
if (ignore_state /= 0) then
write (stdout, '(A,I0,A)') 'Ignoring state ', ignore_state, ' for hopping'
end if
! Determining the initial state
if (irest /= 1) then
call set_current_state(istate_init)
end if
if (nohop == 1) then
write (stdout, '(A,I0)') 'WARNING: nohop=1, Running adiabatic dynamics on state ', istate
end if
deltaE = deltaE / AUtoEV
allocate (nacx(natom, nstate, nstate))
allocate (nacy(natom, nstate, nstate))
allocate (nacz(natom, nstate, nstate))
allocate (nacx_old(natom, nstate, nstate))
allocate (nacy_old(natom, nstate, nstate))
allocate (nacz_old(natom, nstate, nstate))
nacx = 0.0D0
nacy = 0.0D0
nacz = 0.0D0
nacx_old = nacx
nacy_old = nacx
nacz_old = nacx
allocate (en_array(nstate))
allocate (en_array_old(nstate))
en_array = 0.0D0
en_array_old = en_array
! Initialize the history array we use to calculate the Baeck-An couplings
if (inac == 1) then
allocate (en_hist_array(nstate, 4)) !last 3 energies (1: current, 2: n-1, 3: n-2, 4: n-3)
en_hist_array = 0.0D0
allocate (sigma_ba(nstate, nstate)) !this is equivalent to dotproduct, but I will need to store old and new values
sigma_ba = 0.0D0
end if
allocate (tocalc(nstate, nstate))
tocalc = 0
tocalc(istate, istate) = 1
! Compute initial wavefunction and energies,
! used for subsequent determination of TOCALC according to deltaE threshold
! NOTE: We're computing forces as well, since the force_abin() interface
! currently does not allow for computing only energies.
dum_eclas = 0.0D0
dum_fx = 0.0D0; dum_fy = 0.0D0; dum_fz = 0.0D0
call force_clas(dum_fx, dum_fy, dum_fz, x, y, z, dum_eclas, pot)
! open nacme_all.dat if NACME is calculated
if (inac == 0) call nacmefile_init()
! When restarting, initial SH WF was already read from the restart file
if (irest == 0) then
call sh_set_initialwf(istate)
end if
call sh_set_energy_shift(en_array(1))
call sh_select_integrator(integ)
entot0 = en_array(istate) + ekin_v(vx, vy, vz)
entot0 = entot0 * AUtoEV
call set_tocalc()
if (irest == 1) then
call read_nacmrest()
end if
end subroutine sh_init
subroutine read_sh_input(param_unit)
use mod_utils, only: tolower
integer, intent(in) :: param_unit
rewind (param_unit)
read (param_unit, sh)
rewind (param_unit)
integ = tolower(integ)
end subroutine read_sh_input
subroutine print_sh_input()
write (stdout, nml=sh, delim='APOSTROPHE')
end subroutine print_sh_input
subroutine set_current_state(current_state)
integer, intent(in) :: current_state
istate = current_state
end subroutine set_current_state
subroutine check_sh_parameters()
use mod_utils, only: int_positive, int_nonnegative, int_switch, real_nonnegative
use mod_general, only: iknow
use mod_chars, only: chknow
logical :: error
error = .false.
if (integ /= 'euler' .and. integ /= 'rk4' .and. integ /= 'butcher') then
write (stderr, '(A)') 'variable integ must be "euler", "rk4" or "butcher".'
error = .true.
end if
if (integ /= 'butcher') then
write (stderr, '(A)') 'WARNING: variable integ is not "butcher", which is the default and most accurate.'
write (stderr, '(A)') chknow
if (iknow /= 1) error = .true.
end if
if (istate_init > nstate) then
write (stderr, '(A)') 'Initial state > number of computed states.'
error = .true.
end if
! converting input 'couplings' into inac which is used in the code
select case (couplings)
case ('analytic')
inac = 0
write (stdout, '(A)') 'Using analytic ab initio couplings.'
case ('baeck-an')
inac = 1
write (stdout, '(A)') 'Using approximate Baeck-An couplings.'
case ('none')
inac = 2
write (stdout, '(A)') 'Ignoring nonadiabatic couplings.'
case default
write (stderr, '(A)') 'Parameter "couplings" must be "analytic", "baeck-an" or "none".'
error = .true.
end select
! converting input 'velocity_rescaling' into inac which is used in the code
select case (velocity_rescaling)
case ('nac_then_velocity')
adjmom = 0
write (stdout, '(A)') 'Rescaling velocity along the NAC vector after hop.'
write (stdout, '(A)') 'If there is not enough energy, try rescaling along the velocity vector.'
case ('velocity')
adjmom = 1
write (stdout, '(A)') 'Rescaling velocity along the momentum vector after hop.'
case default
write (stderr, '(A)') 'Parameter "velocity_rescaling" must be "nac_then_velocity" or "velocity".'
error = .true.
end select
if (adjmom == 0 .and. inac == 1) then
write (stderr, '(A)') 'Combination of velocity_rescaling="nac_then_velocity" and couplings="baeck-an" is not possible.'
write (stderr, '(A)') 'Velocity cannot be rescaled along NAC when using Baeck-An.'
write (stderr, '(A)') 'Change velocity_rescaling="velocity" to rescale along the velocity vector.'
error = .true.
end if
if (inac == 2 .and. nohop == 0 .and. .not. shwit) then
write (stdout, '(A)') 'WARNING: For simulations without couplings(="none") hopping probability cannot be determined.'
nohop = 1
end if
if (nac_accu1 <= nac_accu2) then
write (stderr, '(A)') 'WARNING: nac_accu1 < nac_accu2'
write (stderr, '(A,I0)') 'Computing NACME only with default accuracy 10^-', nac_accu1
end if
if (decoh_alpha == 0.0D0) then
write (stdout, '(A)') "Turning OFF decoherence correction"
end if
if (ignore_state == istate_init) then
write (stderr, '(A)') 'ERROR: ignore_state == istate_init'
write (stderr, '(A)') 'I cannot start on an ignored state'
error = .true.
end if
if (ignore_state > nstate) then
write (stderr, '(A)') 'Ignored state > number of computed states.'
error = .true.
end if
call int_switch(nohop, 'nohop')
call int_switch(adjmom, 'adjmom')
call int_switch(adjmom, 'revmom')
call int_positive(istate_init, 'istate_init')
call int_positive(nstate, 'nstate')
call int_positive(substep, 'substep')
call int_positive(nac_accu1, 'nac_accu1')
call int_nonnegative(nac_accu2, 'nac_accu2')
call int_nonnegative(ignore_state, 'ignore_state')
call real_nonnegative(decoh_alpha, 'decoh_alpha')
call real_nonnegative(deltaE, 'deltaE')
call real_nonnegative(popthr, 'popthr')
call real_nonnegative(popsumthr, 'popsumthr')
call real_nonnegative(energydifthr, 'energydifthr')
call real_nonnegative(energydriftthr, 'energydriftthr')
call real_nonnegative(dE_S0S1_thr, 'dE_S0S1_thr')
if (error) then
call fatal_error(__FILE__, __LINE__, 'Invalid Surface Hopping parameter(s)')
end if
end subroutine check_sh_parameters
subroutine get_nacm(pot)
character(len=*), intent(in) :: pot
integer :: iost, ist1, ist2
integer :: num_nacme
! In TeraChem SH interface, we already got NACME
if (pot == '_tera_') return
if (pot == '_nai_') return ! for NaI model, the couplings are calcualted together with forces
! Check whether we need to compute any NACME
num_nacme = 0
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
num_nacme = num_nacme + tocalc(ist1, ist2)
end do
end do
if (inac > 0 .or. num_nacme == 0) return
! Calculate NACME using default accuracy
call calc_nacm(pot, nac_accu1)
iost = read_nacm()
if (iost /= 0 .and. nac_accu1 > nac_accu2) then
write (stderr, '(A)') 'WARNING: Some NACME not computed. Trying with decreased accuracy'
write (stderr, '(A,I0)') 'Calling script r.'//trim(pot)//'.nacm with accuracy: 10^-', nac_accu2
call calc_nacm(pot, nac_accu2)
iost = read_nacm()
end if
if (iost /= 0) then
call fatal_error(__FILE__, __LINE__, 'Some NACMEs could not be read')
end if
! we always have to set tocalc because we change it in read_nacm()
call set_tocalc()
end subroutine get_nacm
! In this routine, we decide which gradients and which NACME we need to compute
! TOCALC is a symmetric matrix, upper-triangle defines NACME,
! the diagonal defines gradients
subroutine set_tocalc()
integer :: ist1, ist2
real(DP) :: pop, pop2
tocalc = 0
! The diagonal holds information about gradients that we need
! for SH, we just need the gradient of the current state
tocalc(istate, istate) = 1
! Do not calculate NACME for ADIABATIC dynamics
if (inac == 2) then
return
end if
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
if (abs(en_array(ist1) - en_array(ist2)) < deltaE) then
tocalc(ist1, ist2) = 1
tocalc(ist2, ist1) = 1
end if
end do
end do
if (popthr > 0) then
! COMPUTE NACME only if population of the states is gt.popthr
do ist1 = 1, nstate - 1
pop = get_elpop(ist1)
do ist2 = ist1 + 1, nstate
pop2 = get_elpop(ist2)
if (pop < popthr .and. pop2 < popthr .and. ist1 /= istate .and. ist2 /= istate) then
tocalc(ist1, ist2) = 0
end if
end do
end do
end if
if (ignore_state > 0) then
do ist1 = 1, nstate
tocalc(ist1, ignore_state) = 0
tocalc(ignore_state, ist1) = 0
end do
end if
end subroutine set_tocalc
subroutine write_nacmrest()
use mod_general, only: narchive, it
use mod_qmmm, only: natqm
use mod_utils, only: rename_file, archive_file
integer :: ist1, ist2, iat
integer :: iunit1
character(len=*), parameter :: restart_file = 'restart_sh.bin'
call rename_file(restart_file, trim(restart_file)//'.old')
open (newunit=iunit1, file=restart_file, action='write', status="new", access="sequential", form="unformatted")
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
if (inac == 0) then
do iat = 1, natqm
write (iunit1) nacx(iat, ist1, ist2), &
& nacy(iat, ist1, ist2), &
& nacz(iat, ist1, ist2)
end do
end if
end do
end do
! Printing total energy at t=0, so that we can safely restart
! and we do not break checking for energy drift
! For now, energy jump check is not handled well.
write (iunit1) entot0
if (phase == 1) then
call sh_write_phase_bin(iunit1)
end if
close (iunit1)
if (modulo(it, narchive) == 0) then
call archive_file(restart_file, it)
end if
end subroutine write_nacmrest
subroutine read_nacmrest()
use mod_general, only: it
use mod_qmmm, only: natqm
use mod_utils, only: archive_file
character(len=*), parameter :: restart_file = 'restart_sh.bin'
integer :: iost, ist1, ist2, iat
integer :: iunit1
write (stdout, *) 'Reading Surface Hopping restart data from file '//trim(restart_file)
open (newunit=iunit1, file=restart_file, action="read", status="old", access="sequential", form="unformatted")
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
if (inac == 0) then
do iat = 1, natqm
read (iunit1, iostat=iost) nacx(iat, ist1, ist2), &
& nacy(iat, ist1, ist2), &
& nacz(iat, ist1, ist2)
if (iost /= 0) then
call fatal_error(__FILE__, __LINE__, &
& 'Could not read NACME from restart file '//trim(restart_file))
end if
nacx(iat, ist2, ist1) = -nacx(iat, ist1, ist2)
nacy(iat, ist2, ist1) = -nacy(iat, ist1, ist2)
nacz(iat, ist2, ist1) = -nacz(iat, ist1, ist2)
end do
end if
end do
end do
! Reading total energy at t=0, so that we can safely restart
! and we do not break checking for energy drift
! For now, energy jump check is not handled well.
read (iunit1) entot0
if (phase == 1) then
call sh_read_phase_bin(iunit1)
end if
close (iunit1)
call archive_file(restart_file, it)
end subroutine read_nacmrest
integer function read_nacm()
use mod_qmmm, only: natqm
integer :: iost, ist1, ist2, iat, iunit
iost = 0 ! needed if each tocalc=0
open (newunit=iunit, file='nacm.dat')
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
! TODO: We should have some validation here
! that we're indeed reading NACME between correct states.
if (tocalc(ist1, ist2) == 1) then
do iat = 1, natqm ! reading only for QM atoms
read (iunit, *, IOSTAT=iost) nacx(iat, ist1, ist2), &
& nacy(iat, ist1, ist2), &
& nacz(iat, ist1, ist2)
if (iost == 0) then
tocalc(ist1, ist2) = 0 ! marking as read, useful if we do decreased accuracy
nacx(iat, ist2, ist1) = -nacx(iat, ist1, ist2)
nacy(iat, ist2, ist1) = -nacy(iat, ist1, ist2)
nacz(iat, ist2, ist1) = -nacz(iat, ist1, ist2)
else
close (iunit, status='delete')
write (stderr, '(A,I0,A,I0,A)') 'WARNING: NACME between states ', ist1, ' and ', ist2, ' not read.'
read_nacm = iost
return
end if
end do
! if tocalc
end if
end do
end do
close (iunit, status='delete')
read_nacm = iost
end function read_nacm
subroutine calc_nacm(pot, nac_accu)
use mod_utils, only: toupper
use mod_general, only: it
character(len=*), intent(in) :: pot
integer, intent(in) :: nac_accu
integer :: ist1, ist2, u, itrj
integer :: icmd, istat
character(len=100) :: chsystem
open (newunit=u, file='state.dat')
write (u, '(I2)') nstate
! we print upper triangular part of tocalc matrix to file state.dat
! tocalc(i,j)==1 -> calculate NACME between states i and j
! tocalc(,)==0 -> don't calculate NACME
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
write (u, '(I1,A1)', advance='no') tocalc(ist1, ist2), ' '
end do
end do
close (u)
chsystem = './'//trim(toupper(pot))//'/r.'//trim(pot)//'.nacm '
itrj = 1
write (chsystem, '(A,X,I0,X,I4.3,X,I0,X,A)') trim(chsystem), it, itrj, nac_accu, ' < state.dat'
call execute_command_line(trim(chsystem), exitstat=istat, cmdstat=icmd)
if (icmd /= 0) then
call fatal_error(__FILE__, __LINE__, &
& 'Could not execute command "'//trim(chsystem)//'"')
end if
! We continue on, since we can retry calculation with decreased accuracy
if (istat /= 0) then
write (stderr, *) 'WARNING: Command "'//trim(chsystem)//'" exited with non-zero status'
end if
end subroutine calc_nacm
! Calculate Baeck-An couplings
! Implementation was based on these two articles:
! original by Barbatti: https://doi.org/10.12688/openreseurope.13624.2
! another implementation by Truhlar: https://doi.org/10.1021/acs.jctc.1c01080
! In the numeric implementation, we follow Barbatti and use a higher order formula.
subroutine calc_baeckan(dt)
use mod_general, only: it
integer :: ist1, ist2
real(DP), intent(in) :: dt
real(DP) :: de(4), de2dt2, argument
sigma_ba = 0.0D0
! we don't have sufficient history until step 4
if (it < 4) return
do ist1 = 1, nstate
do ist2 = ist1 + 1, nstate
! If ignore_state is set, we do not calculate sigma (dotproduct) for this state
if (ignore_state == ist1 .or. ignore_state == ist2) cycle
de = en_hist_array(ist2, :) - en_hist_array(ist1, :)
! Second derivative (de2dt2) comes from Eq. 16 from https://doi.org/10.12688/openreseurope.13624.2
de2dt2 = (2.0D0 * de(1) - 5.0D0 * de(2) + 4.0D0 * de(3) - de(4)) / dt**2
argument = de2dt2 / de(1)
if (argument > 0.0D0) then
sigma_ba(ist2, ist1) = dsqrt(argument) / 2.0D0
end if
sigma_ba(ist1, ist2) = -sigma_ba(ist2, ist1)
end do
end do
end subroutine calc_baeckan
! move arrays from new step to old step
subroutine move_vars(vx, vy, vz, vx_old, vy_old, vz_old)
use mod_general, only: natom
real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :)
real(DP), intent(out) :: vx_old(:, :), vy_old(:, :), vz_old(:, :)
integer :: ist1, ist2, iat
do ist1 = 1, nstate
en_array_old(ist1) = en_array(ist1)
if (inac == 0) then ! nedelame, pokud nacitame rovnou dotprodukt
do ist2 = 1, nstate
do iat = 1, natom
nacx_old(iat, ist1, ist2) = nacx(iat, ist1, ist2)
nacy_old(iat, ist1, ist2) = nacy(iat, ist1, ist2)
nacz_old(iat, ist1, ist2) = nacz(iat, ist1, ist2)
end do
end do
end if
end do
! Shift the energy history for Baeck-An couplings
if (inac == 1) then
! Move old energies by 1
en_hist_array(:, 4) = en_hist_array(:, 3)
en_hist_array(:, 3) = en_hist_array(:, 2)
en_hist_array(:, 2) = en_hist_array(:, 1)
! new energy is stored before the couplings are calculated
! I avoided doing the same as with LZSH energy history tracking because then I need to modify force_abin, force_terash and
! every potential in potentials_sh (and all possible future ones). This way it is kept private and does not depend on the
! way energies are calculated.
! See commit: https://github.com/PHOTOX/ABIN/pull/186/commits/918f6837a76ec0f41b575d3ca948253eed2f30cc
end if
vx_old = vx
vy_old = vy
vz_old = vz
end subroutine move_vars
!******************************
! This is the main SH routine !
!******************************
subroutine surfacehop(x, y, z, vx, vy, vz, vx_old, vy_old, vz_old, dt, eclas)
use mod_general, only: natom, nwrite, idebug, it, sim_time, pot
use mod_const, only: AUTOFS
use mod_kinetic, only: ekin_v
real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
real(DP), intent(inout) :: vx(:, :), vy(:, :), vz(:, :)
real(DP), intent(inout) :: vx_old(:, :), vy_old(:, :), vz_old(:, :)
real(DP), intent(inout) :: eclas
real(DP), intent(in) :: dt
! Interpolated quantities
real(DP), dimension(natom, 1) :: vx_int, vy_int, vz_int
real(DP), dimension(natom, 1) :: vx_newint, vy_newint, vz_newint
real(DP), dimension(nstate) :: en_array_int, en_array_newint
real(DP), dimension(natom, nstate, nstate) :: nacx_int, nacy_int, nacz_int
real(DP), dimension(natom, nstate, nstate) :: nacx_newint, nacy_newint, nacz_newint
real(DP), dimension(nstate, nstate) :: dotproduct_int, dotproduct_newint, sigma_ba_old
! Switching probabilities
real(DP) :: t(nstate, nstate)
! Cumulative switching probabilities
real(DP) :: t_tot(nstate, nstate)
real(DP) :: popsum !populations
integer :: ist1, itp
integer :: ihop
! Shortened variable for current state (istate)
! TODO: Rename this!
integer :: ist ! =istate
real(DP) :: fr
real(DP) :: Ekin, dtp
real(DP) :: pop0
! Simulation time in femtoseconds
real(DP) :: stepfs
character(len=500) :: formt
call check_energy(vx_old, vy_old, vz_old, vx, vy, vz)
call check_energydrift(vx, vy, vz)
t_tot = 1.0D0
! First, calculate NACME
if (inac == 0) then ! Analytic ab initio couplings
! For TeraChem MPI / FMS interface, NAC are already computed!
if (pot /= '_tera_' .and. pot /= '_nai_') then
nacx = 0.0D0
nacy = 0.0D0
nacz = 0.0D0
! Compute and read NACME (MOLPRO-SH interface)
call get_nacm(pot)
end if
! TODO: Should we call this with TeraChem?
! I think TC already phases the couplings internally.
call phase_nacme(nacx_old, nacy_old, nacz_old, nacx, nacy, nacz)
else if (inac == 1) then ! Baeck-An couplings
! saving the current energy to the energy history (shifting was already done in previous step in move_vars)
en_hist_array(:, 1) = en_array(:)
sigma_ba_old = sigma_ba ! saving old sigma_ba
call calc_baeckan(dt)
end if
! smaller time step
dtp = dt / substep
! MAIN LOOP
! Smaller time step for electronic population transfer
do itp = 1, substep
ist = istate
! pop0 is later used for Tully's fewest switches
pop0 = get_elpop(ist)
! INTERPOLATION
if ((inac == 0) .or. (inac == 2)) then
fr = real(itp, DP) / real(substep, DP)
call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_newint, vy_newint, vz_newint, &
nacx_newint, nacy_newint, nacz_newint, en_array_newint, &
dotproduct_newint, fr)
fr = real(itp - 1, DP) / real(substep, DP)
call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, &
nacx_int, nacy_int, nacz_int, en_array_int, &
dotproduct_int, fr)
else if (inac == 1) then
fr = real(itp, DP) / real(substep, DP)
call interpolate_ba(vx, vy, vz, vx_old, vy_old, vz_old, vx_newint, vy_newint, vz_newint, &
en_array_newint, dotproduct_newint, sigma_ba, sigma_ba_old, fr)
fr = real(itp - 1, DP) / real(substep, DP)
call interpolate_ba(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, &
en_array_int, dotproduct_int, sigma_ba, sigma_ba_old, fr)
end if
! Integrate electronic wavefunction for one dtp time step
call sh_integrate_wf(en_array_int, en_array_newint, dotproduct_int, dotproduct_newint, dtp)
! Check whether total population is 1.0
popsum = check_popsum()
! Calculate switching probabilities according to Tully's fewest switches
! Probabilities are stored in matrix t
! (although at this point we calculate only the relevant part of the matrix)
call sh_TFS_transmat(dotproduct_int, dotproduct_newint, istate, pop0, t, dtp)
if (idebug > 1) then
stepfs = (sim_time + dtp * itp - substep * dtp) * AUtoFS
call sh_debug_wf(ist, stepfs, t)
end if
do ist1 = 1, nstate
! Cumulative probability over whole big step
! This is only auxiliary variable used for output
t_tot(ist, ist1) = t_tot(ist, ist1) * (1 - t(ist, ist1))
end do
! Should we hop before decoherence?
! Every article says something different.
! Newton-X hops before decoherence.
! HOPPING SECTION
call random_hop(transmat=t, old_state=ist, new_state=ihop)
if (nohop == 1 .and. ihop /= 0) then
write (stdout, '(A,I0)') 'WARNING: Ignoring hop to state ', ihop
end if
! Did HOP occur?
if (nohop /= 1 .and. ihop /= 0) then
! NOTE: Hop can still fail due to insufficient kinetic energy
! ("frustrated hop")
if (adjmom == 0) then
call try_hop_nacme_rescale(vx, vy, vz, ist, ihop, eclas)
else if (adjmom == 1) then
call try_hop_simple_rescale(vx, vy, vz, ist, ihop, eclas)
else
call fatal_error(__FILE__, __LINE__, 'Invalid adjmom value')
end if
! NOTE: We're writing this geometry even for frustrated hops!
! Not sure if that is desired?
call write_hopgeom(x, y, z, old_state=ist, new_state=ihop, timestep=it)
end if
! Apply decoherence correction from Persico et al
if (decoh_alpha > 0) then
Ekin = ekin_v(vx_int, vy_int, vz_int)
if (Ekin > 1.0D-4) then ! Decoherence diverges for zero velocities
call sh_decoherence_correction(en_array_int, decoh_alpha, Ekin, istate, dtp)
end if
end if
!itp loop
end do
if (SHwIT) then
call shwit_check(en_array(1), en_array(2), dE_S0S1_thr, &
vx, vy, vz, eclas)
end if
! set tocalc array for the next step
call set_tocalc()
popsum = check_popsum()
call move_vars(vx, vy, vz, vx_old, vy_old, vz_old)
if (modulo(it, nwrite) == 0) then
call write_sh_output()
end if
contains
subroutine write_sh_output()
use mod_general, only: sim_time
use mod_const, only: ANG, AUTOFS
use mod_files, only: UPOP, UPROB, UPES, UNACME, UDOTPROD
integer :: ist1, ist2, iat
real(DP) :: stepfs
! Simulation time in femtoseconds
stepfs = sim_time * AUtoFS
! Write electronic populations
write (formt, '(A10,I0,A13)') '(F15.2,I3,', nstate, 'F10.5,1F10.7)'
write (UPOP, fmt=formt) stepfs, istate, (get_elpop(ist1), ist1=1, nstate), popsum
! Write hopping probabilities
t_tot = 1 - t_tot ! up to know, t_tot was the probability of not hopping
write (formt, '(A10,I0,A6)') '(F15.2,I3,', nstate, 'F10.5)'
write (UPROB, fmt=formt) stepfs, istate, (t_tot(ist, ist1), ist1=1, nstate)
! Write potential energies to PES.dat
write (formt, '(A7,I0,A7)') '(F15.2,', nstate, 'E20.10)'
write (UPES, fmt=formt) stepfs, (en_array(ist1), ist1=1, nstate)
if (inac == 0) write (UNACME, *) 'Time step:', it
do ist1 = 1, nstate - 1
do ist2 = ist1 + 1, nstate
if (ist1 == 1 .and. ist2 == 2) then
write (UDOTPROD, '(F15.2,E20.10)', advance='no') stepfs, dotproduct_int(ist1, ist2)
else
write (UDOTPROD, '(E20.10)', advance='no') dotproduct_int(ist1, ist2)
end if
if (inac == 0) then
write (UNACME, *) 'NACME between states:', ist1, ist2
do iat = 1, natom
write (UNACME, '(3E20.10)') nacx(iat, ist1, ist2), &
& nacy(iat, ist1, ist2), &
& nacz(iat, ist1, ist2)
end do
end if
end do
end do
write (UDOTPROD, *) ''
end subroutine write_sh_output
end subroutine surfacehop
subroutine phase_nacme(nacx_old, nacy_old, nacz_old, nacx, nacy, nacz)
real(DP), intent(in), dimension(:, :, :) :: nacx_old, nacy_old, nacz_old
real(DP), intent(inout), dimension(:, :, :) :: nacx, nacy, nacz
integer :: ist1, ist2, iat
integer :: natom
real(DP) :: vect_olap
natom = size(nacx, 1)
! Calculating overlap between nacmes
! This is crucial, since NACME vectors can change
! orientation 180 degrees between time steps and we need to correct that
do ist1 = 1, nstate
do ist2 = 1, nstate
vect_olap = 0.0D0
do iat = 1, natom
vect_olap = vect_olap + &
& nacx_old(iat, ist1, ist2) * nacx(iat, ist1, ist2)
vect_olap = vect_olap + &
& nacy_old(iat, ist1, ist2) * nacy(iat, ist1, ist2)
vect_olap = vect_olap + &
& nacz_old(iat, ist1, ist2) * nacz(iat, ist1, ist2)
end do
if (vect_olap < 0) then
do iat = 1, natom
nacx(iat, ist1, ist2) = -nacx(iat, ist1, ist2)
nacy(iat, ist1, ist2) = -nacy(iat, ist1, ist2)
nacz(iat, ist1, ist2) = -nacz(iat, ist1, ist2)
end do
end if
end do
end do
end subroutine phase_nacme
subroutine write_hopgeom(x, y, z, old_state, new_state, timestep)
use mod_system, only: names
use mod_const, only: ANG
real(DP), intent(in), dimension(:, :) :: x, y, z
integer, intent(in) :: old_state, new_state, timestep
character(len=100) :: formt
integer :: natom
integer :: iat, u
natom = size(names)
write (formt, '("hopgeom.",I0,".",I0,".",I0,".xyz")') old_state, new_state, timestep
open (newunit=u, file=trim(formt), action='write')
write (u, *) natom
write (u, *) ''
do iat = 1, natom
write (u, '(A,3ES25.16E3)') names(iat), x(iat, 1) / ANG, y(iat, 1) / ANG, z(iat, 1) / ANG
end do
close (u)
end subroutine write_hopgeom
subroutine random_hop(transmat, old_state, new_state)
use mod_random, only: vranf
! Transition matrix
real(DP), intent(in), dimension(:, :) :: transmat
! Current state index
integer, intent(in) :: old_state
! New state index (0 if no hop occurs)
integer, intent(out) :: new_state
real(DP) :: prob(nstate)
! Pseudorandom number
real(DP) :: rdnum(1)
integer :: ist1
! We return 0 if there's no hop
new_state = 0
! Auxiliary calculations of probabilities on a number line
prob = 0.0D0
! If we are in the ground state, we cannot jump into the ground state :-)
if (old_state /= 1) then
! Probability of jumping from the current state to the ground state
prob(1) = transmat(old_state, 1)
end if
do ist1 = 2, nstate
if (ist1 == old_state) then
prob(ist1) = prob(ist1 - 1)
else
prob(ist1) = prob(ist1 - 1) + transmat(old_state, ist1)
end if
end do
! Get one random number between 0 and 1
call vranf(rdnum, 1)
! Determine, whether we hopped or not
do ist1 = 1, nstate
if (ist1 == old_state) cycle
if (rdnum(1) < prob(ist1)) then
new_state = ist1
exit
end if
end do
end subroutine random_hop
subroutine try_hop_nacme_rescale(vx, vy, vz, instate, outstate, eclas)