forked from alexei-matveev/paragauss-gpl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
integral_calc_quad_2cob3c.f90
4027 lines (3561 loc) · 155 KB
/
integral_calc_quad_2cob3c.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
!
! ParaGauss, a program package for high-performance computations of
! molecular systems
!
! Copyright (C) 2014 T. Belling, T. Grauschopf, S. Krüger,
! F. Nörtemann, M. Staufer, M. Mayer, V. A. Nasluzov, U. Birkenheuer,
! A. Hu, A. V. Matveev, A. V. Shor, M. S. K. Fuchs-Rohr, K. M. Neyman,
! D. I. Ganyushin, T. Kerdcharoen, A. Woiterski, A. B. Gordienko,
! S. Majumder, M. H. i Rotllant, R. Ramakrishnan, G. Dixit,
! A. Nikodem, T. Soini, M. Roderus, N. Rösch
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License version 2 as
! published by the Free Software Foundation [1].
!
! This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! General Public License for more details.
!
! [1] http://www.gnu.org/licenses/gpl-2.0.html
!
! Please see the accompanying LICENSE file for further information.
!
!===================================================================
! Public interface of module
!===================================================================
subroutine integral_calc_quad_2cob3c()
!---------------------------------------------------------------------
!
! Purpose: This routine is the main routine for the
! integral calculation of one quadrupel
! ( unique atom 1, l 1, unique atom 2, l 2 )
! for 2 center orbital and three center integral
! calculation. The corresponding data (including
! the quadrupel to be calculated) are stored in
! int_data_2cob3c_module.
! Contraction and symmetry-adaption are directly
! done within this subroutine.
!
! Subroutine called by: main_slave integral_main_2cob3c
!
! References: Publisher Document: Concepts of Integral Part
!
! Author: TB
! Date: 6/96
!
!===================================================================
! End of public interface of module
!===================================================================
!---------------------------------------------------------------------
! Modifications
!---------------------------------------------------------------------
!
! Modification (Please copy before editing)
! Author: MS
! Date: 3/97
! Description: Necessary changes and supplements for calculation of
! energy gradients
!
! Modification (Please copy before editing)
! Author: MS
! Date: 7/97
! Description: Changes for relativistic and relativistic gradients have
! been added
!
! Modification (Please copy before editing)
! Author: MS
! Date: 8/97
! Description: The symmetryadaption was extended and is now also
! able to treat pseudo 2D irreps
!
! Modification (Please copy before editing)
! Author: Uwe Birkenheuer
! Date: June 1998
! Description: Moving_Unique_Atom concept introduced
! Split_Gradients concept introduced
! Gradients for Model_Density_Approach introduced
!
! Modification (Please copy before editing)
! Author: MM
! Date: 8/97
! Description: Symmetryadaption of spin orbit matrix elements
!
! Modification (Please copy before editing)
! Author: AS
! Date: 11-12/99
! Description: integrals of electrostatic potential are added
!
! Modification (Please copy before editing)
! Modification
! Author: AM
! Date: first: 02.03.1999
! Description: ...
!
! Modification use DLB instead of master/slave for 2cop3c integrals
! Author: AN
! Date: 4/11
! Description: use DLB for scheduling batches of 2cop3c integrals
! remove reporting back of slaves, they get their new
! tasks via DLB
!
! Modification (Please copy before editing)
! Author: ...
! Date: ...
! Description: ...
!---------------------------------------------------------------------
!..............................................................................
! F1 = Sum(n,s) < d/dR psi_n,s | T+V_nuc+V_H - eps_n,s | psi_n,s >
! or
! Sum(n,s) < d/dR psi_n,s | T+V_nuc+V_H+V_X,s - eps_n,s | psi_n,s >
!
! F2 = < rho_tot | d/dR V_nuc >
!
! F3 = < rho_tot | d/dR V_H >
!
! F4 = Sum(s) < rho_s | d/dR V_X,s >
!
!
! Storage
! ~~~~~~~
! options_split_gradients | FALSE FALSE TRUE TRUE
! model_density ! FALSE TRUE FALSE TRUE
! -----------------------------+------------------------------------
! gradient_totalsym(1:n_grads) | F1+F2+F3 F1+F2+F3+F4 F2 F2
! gradient_ob_pulay(1:n_grads) | -- -- F1 F1
! gradient_ch_pulay(1:n_grads) | -- -- F3 F3
! gradient_mda_vxc (1:n_grads) | -- -- -- F4
!..............................................................................
!.......................... STORAGE FOR GRADIENTS .............................
! this a copy of comments from LL_CALCULATE_GRADS:
! << OUTPUT ARRAYS >>
! ===================
! prim_int_2cob_ol_grad ( 1:N) : dsym/dRa < xi_i | 1 | xi_j >
! prim_int_2cob_ol_grad (N+1:M) : dsym/dRb < xi_i | 1 | xi_j >
! << relativistic calculation >>
! prim_int_2cob_kin_grad ( 1:N) : dsym/dRa < xi_i | T | xi_j >
! prim_int_2cob_kin_grad (N+1:M) : dsym/dRb < xi_i | T | xi_j >
! prim_int_2cob_nuc_grad ( ia: ) : dsym/dRa < xi_i | V_nuc | xi_j >
! prim_int_2cob_nuc_grad ( ib: ) : dsym/dRb < xi_i | V_nuc | xi_j >
! prim_int_2cob_nuc_grad ( ic: ) : dsym/dRc < xi_i | V_nuc | xi_j >
! prim_int_2cob_pvsp_grad( ia: ) : dsym/dRa < xi_i | V_pvsp | xi_j >
! prim_int_2cob_pvsp_grad( ib: ) : dsym/dRb < xi_i | V_pvsp | xi_j >
! prim_int_2cob_pvsp_grad( ic: ) : dsym/dRc < xi_i | V_pvsp | xi_j >
! prim_int_3cob_grad ( ia: ) : dsym/dRa < xi_i | V_H | xi_j >
! prim_int_3cob_grad ( ib: ) : dsym/dRb < xi_i | V_H | xi_j >
! prim_int_3cob_grad ( ic: ) : dsym/dRc < xi_i | V_H | xi_j >
! << non-relativistic calculation with total gradients >>
! prim_int_3cob_grad ( ia: ) : dsym/dRa < xi_i | T + V_nuc + V_H | xi_j >
! prim_int_3cob_grad ( ib: ) : dsym/dRb < xi_i | T + V_nuc + V_H | xi_j >
! prim_int_3cob_grad ( ic: ) : dsym/dRc < xi_i | V_nuc + V_H | xi_j >
! << non-relativistic calculation with split gradients >>
! prim_int_2cob_ks_grad ( 1:N) : dsym/dRa < xi_i | T + V_nuc + V_H | xi_j >
! prim_int_2cob_ks_grad (N+1:M) : dsym/dRb < xi_i | T + V_nuc + V_H | xi_j >
! prim_int_3cob_nuc_grad ( ic: ) : dsym/dRc < xi_i | V_nuc | xi_j >
! prim_int_3cob_grad ( ic: ) : dsym/dRc < xi_i | V_H | xi_j >
!------------ Modules used ------------------------------------
! use CPU_TIME:
! define FPP_TIMERS 2
# include "def.h"
! define _NUSD_BY_CALC_3CENTER
use type_module ! type specification parameters
use datatype
use strings, only: itoa
use timer_module
use time_module
use integralpar_module, only: integralpar_gradients &
, integralpar_dervs &
, integralpar_pot_for_secderiv &
, integralpar_totalsymmetric &
, integralpar_3cob_grad &
, integralpar_3c_co &
, integralpar_3c_co_resp &
, integralpar_3c_r2_pvsp &
, integralpar_3c_r2_pvxp &
, integralpar_3c_xc &
, integralpar_3c_rcoul_pvsp &
, integralpar_3c_rcoul_pvxp &
, integralpar_send_3c &
, integralpar_i_int_part &
, integralpar_relativistic &
, integralpar_2cob_pvec &
, integralpar_2cob_potential &
, integralpar_2cob_field &
, integralpar_solv_grad &
, integralpar_Q_solv_grad &
, integralpar_2cob_pc_grad &
, integralpar_2cob_X_grad &
, integralpar_2cob_ipd_grad &
, integralpar_2cob_kin &
, integralpar_2cob_kin_grad &
, integralpar_2cob_nuc &
, integralpar_2cob_nuc_grad &
, integralpar_2cob_pvsp &
, integralpar_2cob_pvsp_grad &
, integralpar_2cob_pvxp &
, integralpar_2cob_ol &
, integralpar_2cob_ol_grad &
, integralpar_cpks_contribs &
, integralpar_cpksdervs &
, integralpar_pseudo
#ifdef WITH_EFP
use integralpar_module, only: integralpar_2cob_pc &
, integralpar_2cob_X &
, integralpar_2cob_field_efp &
, integralpar_efp_gradients
#endif
use int_data_2cob3c_module
use integral_calc_quad_module, only: &
integral_calc_quad_densmat, quad_pmat, quad_wmat &
, integral_calc_quad_close &
, IDONT_CONTRACT, IDONT_SUM_OVER_PARTNERS &
, contrsym2c, contrsym3c, contr3c &
, quad_density_mat
use output_module
use int_send_2cob3c_module, only : int_send_2cob3c_send
use int_send_2cob3c_spor_module, only : int_send_2cob3c_spor_send
use iounitadmin_module
use gradient_data_module, only: gradient_totalsym &
, gradient_data_n_gradients &
, gradient_data_n_spin_gradients &
, dervs_totalsym &
, gradient_ob_pulay &
, gradient_ch_pulay &
, gradient_mda_vxc &
, gradient_index
#ifdef WITH_EPE
use gradient_data_module, only: calc_cluster_epe_energy, cluster_epe_energy
#endif
use unique_atom_module, only: unique_atom_symequivatoms_type &
, unique_atom_partner_type &
, unique_atom_sa_int_type &
, unique_atom_symequiv &
, unique_atoms &
, N_unique_atoms &
, pseudopot_present
use symmetry_data_module
use comm_module
use msgtag_module
use options_module, only: options_split_gradients, options_n_spin, &
xcmode_model_density, xcmode_extended_mda, &
options_xcmode, options_spin_orbit
#ifdef WITH_EPE
use ewaldpc_module,only: ewpc_n
#endif
use pointcharge_module, only: totsym_grad_pc_length
use point_dqo_module, only: totsym_grad_dip_length, &
totsym_grad_quad_length,totsym_grad_oct_length, totsym_grad_rep_length
use induced_dipoles_module, only: totsym_grad_idip_length
use elec_static_field_module, only: totsym_field_length
use solv_cavity_module, only: with_pc, fixed_pc
use spin_orbit_module, only: is_on,op_FitTrafo
use error_module, only: MyID !<<< debug
use calc3c_switches!!!,only: integralpar_cpksdervs,integralpar_2dervs,integralpar_dervs
USE_MEMLOG
use efield_module, only: efield_applied
#ifdef WITH_RESPONSE
use int_send_3c_resp, only: int_send_3c_resp_save
#endif
use interfaces ! FIXME: need it here?
#ifdef WITH_DFTPU
use dft_plus_u_module, only: dft_plus_u_grad, dft_plus_u_in_use &! for DFT+U gradients
, dft_plus_u_mo_in_use
#endif
implicit none
!
! interfaces moved to
! modules/interfaces.f90
!
!------------ Declaration of local variables ------------------
type(unique_atom_symequivatoms_type), pointer :: symequivatoms
! different sym.equiv. atom pairs to be processed
integer(kind=i4_kind) :: i_symequiv
! index of not-symmetry-equivalent connection vector
real(r8_kind) :: weight
! weight of not-symmetry-equivalent connection vector
integer(kind=i4_kind) :: i_ua1, i_ua2, i_l1, i_l2
! local vars for for quadrupel%ua1, quadrupel%ua2, ...
integer(kind=i4_kind) :: i_ea1, i_ea2
! loop indices for equal atoms
logical :: split_gradients, model_density
real(r8_kind), parameter :: zero = 0.0_r8_kind
logical :: &
diagonal_ua, & ! if ( ua1 == ua2 )
diagonal_ea, & ! if ( ua1 == ua2 && ea1 = ea2 )
diagonal_sh ! if ( ua1 == ua2 && ea1 = ea2 && l1 == l2 )
logical :: do_sum_up
#ifdef WITH_EFP
logical :: do_sum_up_efp
#endif
#ifdef FPP_TIMERS
real(r8_kind), parameter :: interval=300.0
real(r8_kind) :: last_time=-2*interval
FPP_TIMER_DECL(i2c3c)
FPP_TIMER_DECL(oli)
FPP_TIMER_DECL(olg)
FPP_TIMER_DECL(shgi)
FPP_TIMER_DECL(cntr3)
FPP_TIMER_DECL(sym3c)
FPP_TIMER_DECL(sumup)
FPP_TIMER_DECL(tomo)
#endif
!--------------------------------------------------------------
!------------ Executable code ---------------------------------
FPP_TIMER_START(i2c3c)
i_ua1 = quadrupel%ua1
i_l1 = quadrupel%l1
i_ua2 = quadrupel%ua2
i_l2 = quadrupel%l2
call dtrace("start with quadrupel "//itoa(i_ua1)//" " &
//itoa(i_l1) //" " &
//itoa(i_ua2)//" " &
//itoa(i_l2) &
)
call dtrace("setup")
call setup()
split_gradients = options_split_gradients()
model_density = options_xcmode() == xcmode_model_density .or. &
options_xcmode() == xcmode_extended_mda
! decide whether or not this is a "property" run that will
! require taking traces with density matrix (see sum_up_gradient)
do_sum_up = integralpar_gradients .or. integralpar_pot_for_secderiv
#ifdef WITH_EFP
do_sum_up = do_sum_up .or. integralpar_2cob_field_efp .or. &
integralpar_2cob_pc .or. integralpar_2cob_X .or. integralpar_efp_gradients
do_sum_up_efp = integralpar_2cob_field_efp.or. &
integralpar_2cob_pc .or. integralpar_2cob_X .or. integralpar_efp_gradients
#endif
if( do_sum_up )then
! if this is a property run, then prepare the section of the
! densmat (and energy-weighted densmat) corresponding to the current
! quadrupel:
call integral_calc_quad_densmat(i_ua1, i_l1, i_ua2, i_l2, quad_pmat, quad_wmat)
endif
! allocate primitives now here, since same for all ea's
! call allocate_primitives()
! MOVED TO int_data_2cob3c_setup()
!
! FIXME: do the same with contracted
! currently cont_* are allocated in
! contract_2center() and contract_3center()
! for each ea
totalsymmetric: if (integralpar_totalsymmetric) then
! In case only total symmertic integrals are calculated,
! integrals are only necessary between selected equal atoms
symequivatoms => unique_atom_symequiv(quadrupel%ua1,quadrupel%ua2)
! loop over not symmetry equivalent pairs of equal atoms
symequiv: do i_symequiv = 1, symequivatoms%n_equiv_dist
call dtrace("processing pair of non-symmetryequivalent atoms: "//itoa(i_symequiv))
i_ea1 = 1
i_ea2 = symequivatoms%index_eq_atom2(i_symequiv)
center1 = ua1%position(:,i_ea1)
center2 = ua2%position(:,i_ea2)
! weight of a pair (number of symm-equiv distances):
weight = symequivatoms%weight(i_symequiv)
diagonal_ua = ( i_ua1 == i_ua2 )
diagonal_ea = ( diagonal_ua .and. (i_ea1 == i_ea2) )
diagonal_sh = ( diagonal_ea .and. (i_l1 == i_l2 ) )
call dtrace("index of secound equal atom: "//itoa(i_ea2) )
! allocation used to have been done inside of
! call calc_primitives_and_contract()
! it was moved outside the loop over ea-ea distances
! some ?buggy? ??_calculate_*() progs reguire re-initialization
! of primitives to zero
! FIXME: better fix the actual subs not to depend on input state
call allocate_primitives('initialize')
! calculate primitive integrals and contract them in non_relativistic
! case or else map contracted integrals to primitives
call dtrace("calc_primitives")
#ifdef no_cpks_coul_grads
if(allocated(prim_int_dens_mat))then
! it is allocated in the (first?) gradient run
! for computing (explicit?) gradients due to fit-functions
! See int_data_2cob3c_module.
! only in first gradient run of the sec-der calculations:
! Compute the section of the density matrix in (uncontracted?) basis
! that corresponds to the curent quadrupel.
call quad_density_mat(i_ua1,i_ea1,i_l1,i_ua2,i_ea2,i_l2,weight,prim_int_dens_mat)
endif
#endif
call calc_primitives(i_ua1,i_ea1,i_l1,i_ua2,i_ea2,i_l2)
! add contribution of one not symmetry equivalent pair
! of equal atom to symmetry adapted integrals
call dtrace("adding results to symmetry adapted arrays(a)")
call contract_and_symadapt(i_ua1,i_ea1,i_l1,i_ua2,i_ea2,i_l2,0) ! sum over partners totalsym
call dtrace("1 pair of non-symmetryequivalent atoms done ")
enddo symequiv
else ! now not total symmetric
eqal_atom_1: do i_ea1 = 1, ua1%N_equal_atoms
eqal_atom_2: do i_ea2 = 1, ua2%N_equal_atoms
call dtrace("processing pair of equal atoms "//itoa(i_ea1)//" and "//itoa(i_ea2))
center1 = ua1%position(:,i_ea1)
center2 = ua2%position(:,i_ea2)
! weight of a pair (all ea1-ea2 distances):
weight = 1.0_r8_kind
call allocate_primitives('initialize')
! calculate primitive integrals and contract them in non_relativistic
! case or else map contracted integrals to primitives
call dtrace("calc_primitives_and_contract")
call calc_primitives(i_ua1,i_ea1,i_l1,i_ua2,i_ea2,i_l2)
! add contribution of one pair
! of equal atom to symmetry adapted integrals
call dtrace("adding results to symmetry adapted arrays(b)")
! not a totally symmetric case:
call contract_and_symadapt(i_ua1, i_ea1, i_l1, i_ua2, i_ea2, i_l2, IDONT_SUM_OVER_PARTNERS)
#ifdef WITH_RESPONSE
if (integralpar_3c_co_resp) then
call symadapt_add_nottotalsym_co(i_l1, i_l2, i_ea1, i_ea2, contr3c(prim_int_3c_co))
end if
#endif
call dtrace("1 pair of equal atoms done")
end do eqal_atom_2
end do eqal_atom_1
endif totalsymmetric !/else
! FIXME: also in not-totalsymmetric case?:
if(.not.options_spin_orbit)then
! in the case of groups with 2D pseudo irreps we have to make
! an adittional symmetry adaption
call symadp_pseudo2D() !!! the only call
endif
call dtrace("calculation done")
if(do_sum_up) then
#ifdef WITH_EFP
if(.not. do_sum_up_efp) then
#endif
call dtrace("sum_up_gradient")
call sum_up_gradient()
#ifdef WITH_DFTPU
!*************************************************************************************
! DFT+U Gradient subroutine is called here
!*************************************************************************************
! DFT+U Gradient is calculated only when both "integralpar_3cob_grad" and
! "dft_plus_u_in_use" are "t". When only "dft_plus_u_in_use" was used, gradient
! calculations with both DFT+U and SOLVATION terminated with error.
if( integralpar_3cob_grad .and. (dft_plus_u_in_use .or. dft_plus_u_mo_in_use) )then
call dft_plus_u_grad(i_ua1, i_l1, i_ua2, i_l2, symadapt_int_2cob_ol_grad, gradient_totalsym)
end if
#endif
#ifdef WITH_EFP
else
call dtrace("sum_up_grads_on_efp")
call sum_up_grads_on_efp()
end if
#endif
! free the storage in integral_calc_quad_module
! that was allocated by integral_calc_quad_densmat(..)
! before proceeding to the next quadrupel:
call integral_calc_quad_close()
endif
if ( integralpar_send_3c .and. .not. integralpar_pot_for_secderiv) then
call start_timer(timer_int_send_2cob3c(integralpar_i_int_part))
call dtrace("int_send_2cob3c_send")
if( .not. options_spin_orbit )then
call int_send_2cob3c_send()
else
call int_send_2cob3c_spor_send()
endif
call stop_timer(timer_int_send_2cob3c(integralpar_i_int_part))
endif
! store relativistic gradients/derivatives for future transformations:
if( integralpar_relativistic .and..not. options_spin_orbit )then
call relstore(i_ua1,i_l1,i_ua2,i_l2)
endif
#ifdef WITH_RESPONSE
if ( integralpar_3c_co_resp ) then
call start_timer(timer_int_send_2cob3c(integralpar_i_int_part))
call dtrace("int_send_2cob3c_send")
call int_send_3c_resp_save(quadrupel%ua1, quadrupel%ua2, &
quadrupel%l1, quadrupel%l2, symadapt_int_3c_co_resp)
call stop_timer(timer_int_send_2cob3c(integralpar_i_int_part))
endif
#endif
call dtrace("shutdown")
call shutdown()
call dtrace("end")
FPP_TIMER_STOP(i2c3c)
#ifdef FPP_TIMERS
if( FPP_TIMER_VALUE(i2c3c) - last_time > interval )then
last_time = FPP_TIMER_VALUE(i2c3c)
print*,MyID,'TIMER: i2c3c tot =',FPP_TIMER_VALUE(i2c3c)
print*,MyID,'TIMER: old ints =',FPP_TIMER_VALUE(oli)
if( integralpar_gradients )then
print*,MyID,'TIMER: old grads=',FPP_TIMER_VALUE(olg)
endif
print*,MyID,'TIMER: shgi ints =',FPP_TIMER_VALUE(shgi)
print*,MyID,'TIMER: contr 3c =',FPP_TIMER_VALUE(cntr3)
print*,MyID,'TIMER: symadp 3c =',FPP_TIMER_VALUE(sym3c)
if( integralpar_gradients )then
print*,MyID,'TIMER: sum up =',FPP_TIMER_VALUE(sumup)
print*,MyID,'TIMER: |- to MO =',FPP_TIMER_VALUE(tomo)
endif
endif
#endif
!--------------------------------------------------------------
!------------ Private Subroutines -----------------------------
contains ! integral_calc_quad_2cob3c
subroutine dtrace(msg)
implicit none
character(len=*), intent(in) :: msg
! *** end of interface ***
if ( output_int_loops ) &
call write_to_output_units(MyID//"integral_calc_quad_2cob3c: "//msg)
end subroutine dtrace
!*************************************************************
subroutine setup
! Purpose: setup routines of various modules
!------------ Executable code ------------------------------
call stop_timer(timer_int_idle_2cob3c(integralpar_i_int_part))
call init_timer(timer_int_calc_2cob3c(integralpar_i_int_part))
call start_timer(timer_int_calc_2cob3c(integralpar_i_int_part))
call start_timer(timer_int_quadrupel_2cob3c(integralpar_i_int_part))
call int_data_2cob3c_setup()
end subroutine setup
subroutine shutdown
! Purpose: shutdown routines of various modules
!------------ Executable code ------------------------------
call int_data_2cob3c_shutdown()
call stop_timer(timer_int_calc_2cob3c(integralpar_i_int_part))
call stop_timer(timer_int_quadrupel_2cob3c(integralpar_i_int_part))
call start_timer(timer_int_idle_2cob3c(integralpar_i_int_part))
call timer_small_to_large( &
timer_int_calc_2cob3c(integralpar_i_int_part), &
timer_int_calcsum_2cob3c(integralpar_i_int_part) )
end subroutine shutdown
!**************************************************************
!**************************************************************
!
! allocate_primitives()
! deallocate_primitives()
! deallocate_contracted()
! all moved to int_data_2cob3c_module
!
!**************************************************************
!**************************************************************
! These were moved (with slight modifications) to
! integral_calc_quad_module:
!
! subroutine contrsym2c *
! subroutine contrsym3c *
! subroutine contract_2center
! subroutine uncontract_2center
! subroutine contract_3center
! function contr2c *
! function contr3c
!
! Only those marked by star are exported back to this sub.
!**************************************************************
!**************************************************************
subroutine calc_primitives(U1,E1,L1,U2,E2,L2)
! calculate primitive integrals and do contraction
use dip_prim_module, only: ll_momnt
use spin_orbit_module !<<< options only
use cpksdervs_matrices
use calc3c_switches
use interfaces
use shgi_cntrl, only: IPSEU
use shgi_cntrl, only: INUSD
implicit none
integer(i4_kind), intent(in) :: U1,E1,L1,U2,E2,L2
! uniqe atom, equiv. atom, and ang momenta for center 1 and 2
!------------ Declaration of local variables ---------------
integer(kind=i4_kind), parameter :: SS=0,LS=1,LL=2
integer(kind=i4_kind) :: ILL ! one of the above
integer(kind=i8_kind) :: IMODE
!------------ Executable code ------------------------------
ILL = LS
if( L1==0 .and. L2==0 ) ILL = SS
if( L1/=0 .and. L2/=0 ) ILL = LL
IMODE = 0
! no more supported:
!ifdef WITH_OLD_PSEUDO
! WARN('WITH_OLD_PSEUDO')
! if( pseudopot_present ) IMODE = IOR(IMODE,IPSEU)
! ! otherwise it will be computed by SHGI
!endif
#ifdef _NUSD_BY_CALC_3CENTER
! tell calc_3center() explicitly to add sec der of nuc:
if( integralpar_dervs ) IMODE = IOR(IMODE,INUSD)
#endif
! calculate primitive integrals
call start_timer(timer_int_prim_2cob3c(integralpar_i_int_part))
#ifdef WITH_EFP
if (.not.integralpar_gradients .and. .not.integralpar_efp_gradients) then
#else
if (.not.integralpar_gradients) then
#endif
FPP_TIMER_START(oli)
#ifdef NEW_INTEGRALS
call start_timer(timer_new_ll_integral)
call calc_3center(U1, L1, U2, L2, IMODE=IMODE)
call stop_timer(timer_new_ll_integral)
#endif
! FIXME: is this debugging?
if(cpks_noco) then
prim_int_3c_co = 0.0_r8_kind
endif
FPP_TIMER_STOP(oli)
if (integralpar_2cob_pvec) then
call ll_momnt(&
& ua2_basis%exponents, ua1_basis%exponents,& !from int_data_2cob_
& L2,L1,& !...
& center2, center1, & !...
& prim_2cv_ints(int_prim_pvec)%prim &
& )
endif
else ! ========= GRADIENT PART ==========
FPP_TIMER_START(olg)
select case(ILL)
case (SS)
call dtrace("calc_primitives: calling ss_calculate_grads")
#ifndef NEW_INTEGRALS
call start_timer(timer_ss_integral)
call ss_calculate_grads(i_ea1,i_ea2,IMODE)
call stop_timer(timer_ss_integral)
#else
call start_timer(timer_new_ss_integral)
FPP_TIMER_START(t_calc_3center)
#if 1 /* commented to see memory use */
call calc_3center(U1,L1,U2,L2, &
equalb=i_ea2,equala=i_ea1,IMODE=IMODE)
#endif
FPP_TIMER_STOP(t_calc_3center)
call stop_timer(timer_new_ss_integral)
#endif
case(LS)
call dtrace("calc_primitives: calling ls_calculate_grads")
#ifndef NEW_INTEGRALS
call start_timer(timer_ls_integral)
call ls_calculate_grads(U1,U2,i_ea2,L1,L2,IMODE)
call stop_timer(timer_ls_integral)
#else
call start_timer(timer_new_ls_integral)
FPP_TIMER_START(t_calc_3center)
#if 1 /* commented to see memory use */
call calc_3center(U1,L1,U2,L2,equalb=i_ea2,IMODE=IMODE)
#endif
FPP_TIMER_STOP(t_calc_3center)
call stop_timer(timer_new_ls_integral)
#endif
case (LL)
call dtrace("calc_primitives: calling ll_calculate_grads")
#ifndef NEW_INTEGRALS
call start_timer(timer_ll_integral)
call ll_calculate_grads(U1,U2,i_ea2,L1,L2,IMODE)
call stop_timer(timer_ll_integral)
#else
call start_timer(timer_new_ll_integral)
FPP_TIMER_START(t_calc_3center)
#if 1 /* commented to see memory use */
call calc_3center(U1,L1,U2,L2,equalb=i_ea2,IMODE=IMODE)
#endif
FPP_TIMER_STOP(t_calc_3center)
call stop_timer(timer_new_ll_integral)
#endif
case default
call error_handler( &
"calc_quad_2cob3c: calc_primitives:&
& a wonder")
end select
FPP_TIMER_STOP(olg)
endif
DPRINT 'call calc_prim_shgi(',U1,E1,L1,U2,E2,L2,')'
FPP_TIMER_START(shgi)
if(integralpar_2cob_potential.or.integralpar_2cob_field) then
call calc_prim_pot_and_field(U1,E1,L1,U2,E2,L2)
else if (integralpar_solv_grad .or.integralpar_Q_solv_grad)then
call calc_prim_solv(U1,E1,L1,U2,E2,L2)
else
#ifdef WITH_EFP
if(.not. integralpar_2cob_pc .and. &
.not. integralpar_2cob_X .and. &
.not.integralpar_efp_gradients)then
#endif
call calc_prim_shgi(U1,E1,L1,U2,E2,L2)
#ifdef WITH_EFP
endif
#endif
if( &
#ifdef WITH_EFP
integralpar_2cob_pc .or. &
integralpar_2cob_X .or. &
#endif
integralpar_2cob_pc_grad.or. &
integralpar_2cob_X_grad .or. &
integralpar_2cob_ipd_grad) then
call dtrace("calc_primitives: calling calc_prim_shgi_X")
call calc_prim_shgi_X(U1,E1,L1,U2,E2,L2)
end if
! calculation of forces in presence of electric field
if (efield_applied() .and. integralpar_gradients) then
call calc_grad_efield(U1,E1,L1,U2,E2,L2)
endif
endif
FPP_TIMER_STOP(shgi)
call dtrace("calc_primitives: primitive integrals done")
call stop_timer(timer_int_prim_2cob3c(integralpar_i_int_part))
end subroutine calc_primitives
!**************************************************************
subroutine contract_and_symadapt(U1, E1, L1, U2, E2, L2, imode)
! calculate primitive integrals and do contraction
use fit_coeff_module, only: fit_coeff_n_ch
implicit none
integer(i4_kind), intent(in) :: U1,E1,L1,U2,E2,L2
integer(i4_kind), intent(in) :: imode
! uniqe atom, equiv. atom, and ang momenta for center 1 and 2
!------------ Declaration of local variables ---------------
integer(kind=i4_kind) :: GRAD_DIM,SPIN_GRAD_DIM
integer(kind=i4_kind) :: i_grad,k2dr
integer(kind=i4_kind) :: irel2c
integer(kind=i4_kind) :: irel3c
integer(kind=i4_kind) :: n_ff_ch ! fit_coeff_n_ch()
!------------ Executable code ------------------------------
call start_timer(timer_int_cont_2cob3c(integralpar_i_int_part))
if (options_spin_orbit) then
! contraction of 3c_co is done inside:
call symadapt_add_totalsym_spor(weight)
GOTO 999 ! finalize and exit
endif
! === NORMAL INTEGRALS ===
! default is imode, in relativistic calculations
! some 2c/3c integrals require special treatment:
irel2c = imode
if ( integralpar_relativistic ) irel2c = IOR(irel2c,IDONT_CONTRACT)
irel3c = imode
if ( is_on(op_FitTrafo) ) irel3c = IOR(irel3c,IDONT_CONTRACT)
if ( integralpar_2cob_kin ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_kin &
,symadapt_int_2cob_kin &
,irel2c &
)
endif
if ( integralpar_2cob_nuc ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_nuc &
,symadapt_int_2cob_nuc &
,irel2c &
)
if((pseudopot_present.and.integralpar_relativistic) .or. integralpar_pseudo) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_nuc_pseudo &
,symadapt_int_2cob_nuc_pseudo &
,irel2c &
)
endif
endif
if ( integralpar_2cob_pvsp ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_pvsp &
,symadapt_int_2cob_pvsp &
,irel2c &
)
endif
if ( integralpar_2cob_ol ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_ol &
,symadapt_int_2cob_ol &
,irel2c &
)
endif
if ( integralpar_3c_xc ) then
call contrsym3c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_3c_xc &
,symadapt_int_3c_xc &
,irel3c &
)
endif
if ( integralpar_3c_co ) then
!
! FIXME: int TDDFT (response) calculations, the shape of
! the primitive fit integrals, prim_int_3c_co, deviates
! from normal. The middle axis has the length of a total
! number of fit functions including the non-totally symmetric
! once. As the totally symmetric irrep is usually the first
! the leading once are those used for density expansion:
!
n_ff_ch = fit_coeff_n_ch()
call contrsym3c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_3c_co(:, :, :n_ff_ch, :, :) &
,symadapt_int_3c_co &
,irel3c &
)
endif
if ( integralpar_2cob_potential ) then
call contrsym3c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_poten &
,symadapt_int_2cob_poten &
,imode &
)
endif
if ( integralpar_2cob_field ) then
call contrsym3c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_field &
,symadapt_int_2cob_field &
,imode &
)
endif
! ========================
! === GRADIENTS ===
grad_dim=size(symadapt_int_2cob_ol_grad,1)
! split_gradients = options_split_gradients()
! model_density = options_xcmode() == xcmode_model_density .or. &
! options_xcmode() == xcmode_extended_mda
if (model_density) then
spin_grad_dim = grad_dim * options_n_spin()
else
spin_grad_dim = grad_dim
endif
#ifdef WITH_EPE
if ( integralpar_3cob_grad .and. ewpc_n.ne.0.and.calc_cluster_epe_energy) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_3cob_epe%m &
,symadapt_int_3cob_epe(1,:) &
,imode &
)
endif
#endif
do i_grad=1,GRAD_DIM
if ( integralpar_2cob_ol_grad ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_ol_grad(i_grad)%m &
,symadapt_int_2cob_ol_grad(i_grad,:) &
,imode &
)
if( integralpar_relativistic ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_ol_grad(i_grad)%m &
,symadapt_int_2cob_ol_rel_grad(i_grad,:) &
,irel2c &
)
endif
endif
if ( integralpar_2cob_kin_grad ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_kin_grad(i_grad)%m &
,symadapt_int_2cob_kin_grad(i_grad,:) &
,irel2c &
)
endif
enddo
do i_grad=1,gradient_data_n_gradients
if ( integralpar_2cob_nuc_grad ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_nuc_grad(i_grad)%m &
,symadapt_int_2cob_nuc_grad(i_grad,:) &
,irel2c &
)
if((pseudopot_present.and.integralpar_relativistic) .or. integralpar_pseudo) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_pseudo_grad(i_grad)%m &
,symadapt_int_2cob_pseudo_grad(i_grad,:) &
,irel2c &
)
endif
endif
if ( integralpar_2cob_pvsp_grad ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_2cob_pvsp_grad(i_grad)%m &
,symadapt_int_2cob_pvsp_grad(i_grad,:) &
,irel2c &
)
endif
if ( integralpar_solv_grad .and. integralpar_cpksdervs ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_3cob_solv_grad(i_grad)%m &
,symadapt_int_3cob_solv_grad(i_grad,:) &
,imode &
)
endif
if ( integralpar_3cob_grad ) then
call contrsym2c( &
U1,E1,L1,U2,E2,L2,weight &
,prim_int_3cob_grad(i_grad)%m &