9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-20 17:53:25 +02:00

Merge pull request #32 from QuantumPackage/dev-stable

Dev stable
This commit is contained in:
AbdAmmar 2024-02-19 10:47:41 +01:00 committed by GitHub
commit 6b9649fc2c
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
28 changed files with 2127 additions and 46 deletions

View File

@ -1,3 +1,7 @@
**Important**: The Intel ifx compiler is not able to produce correct
executables for Quantum Package. Please use ifort as long as you can, and
consider switching to gfortran in the long term.
# Quantum Package 2.2
<!--- img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250" --->

62
config/gfortran_mkl.cfg Normal file
View File

@ -0,0 +1,62 @@
# Common flags
##############
#
# -ffree-line-length-none : Needed for IRPF90 which produces long lines
# -lblas -llapack : Link with libblas and liblapack libraries provided by the system
# -I . : Include the curent directory (Mandatory)
#
# --ninja : Allow the utilisation of ninja. (Mandatory)
# --align=32 : Align all provided arrays on a 32-byte boundary
#
#
[COMMON]
FC : gfortran -ffree-line-length-none -I . -mavx -g -fPIC -std=legacy
LAPACK_LIB : -I${MKLROOT}/include -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_core -lpthread -lm -ldl -lmkl_gnu_thread -lgomp -fopenmp
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED
# Global options
################
#
# 1 : Activate
# 0 : Deactivate
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags
####################
#
# -Ofast : Disregard strict standards compliance. Enables all -O3 optimizations.
# It also enables optimizations that are not valid
# for all standard-compliant programs. It turns on
# -ffast-math and the Fortran-specific
# -fno-protect-parens and -fstack-arrays.
[OPT]
FCFLAGS : -Ofast -mavx
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast
# Debugging flags
#################
#
# -fcheck=all : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
#
[DEBUG]
FCFLAGS : -fcheck=all -g
# OpenMP flags
#################
#
[OPENMP]
FC : -fopenmp
IRPF90_FLAGS : --openmp

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
FC : ifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : mpiifort -fpic
FC : mpiifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
FC : ifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --define=WITHOUT_TRAILZ --define=WITHOUT_SHIFTRL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
FC : ifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : mpiifort -fpic
FC : mpiifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
FC : ifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
FC : ifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : mpiifort -fpic
FC : mpiifort -fpic -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL

View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic -diag-disable 5462
FC : ifort -fpic -diag-disable=5462 -diag-disable=10448
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=64 -DINTEL

View File

@ -58,17 +58,32 @@ let int_of_atom_id : atom_id -> int = fun x -> x
let float_of_distance : float StringMap.t -> distance -> float =
fun map -> function
| Value x -> x
| Label s -> StringMap.find s map
| Label s -> begin
try StringMap.find s map with
| Not_found ->
Printf.sprintf "Zmatrix error: distance %s undefined" s
|> failwith
end
let float_of_angle : float StringMap.t -> angle -> float =
fun map -> function
| Value x -> x
| Label s -> StringMap.find s map
| Label s -> begin
try StringMap.find s map with
| Not_found ->
Printf.sprintf "Zmatrix error: angle %s undefined" s
|> failwith
end
let float_of_dihedral : float StringMap.t -> dihedral -> float =
fun map -> function
| Value x -> x
| Label s -> StringMap.find s map
| Label s -> begin
try StringMap.find s map with
| Not_found ->
Printf.sprintf "Zmatrix error: dihedral %s undefined" s
|> failwith
end
type line =

View File

@ -224,7 +224,7 @@
subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
BEGIN_DOC
! Computes the following integral :
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ]
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x^2 ]
! needed for the dipole and those things
END_DOC
implicit none

View File

@ -4,13 +4,15 @@ casscf
|CASSCF| program with the CIPSI algorithm.
Example of inputs
-----------------
Example of inputs for GROUND STATE calculations
-----------------------------------------------
NOTICE :: FOR EXCITED STATES CALCULATIONS SEE THE FILE "example_casscf_multistate.sh"
a) Small active space : standard CASSCF
---------------------------------------
Let's do O2 (triplet) in aug-cc-pvdz with the following geometry (xyz format, Bohr units)
3
2
O 0.0000000000 0.0000000000 -1.1408000000
O 0.0000000000 0.0000000000 1.1408000000
@ -45,3 +47,4 @@ qp set casscf_cipsi small_active_space False
qp run casscf | tee ${EZFIO_FILE}.casscf_large.out
# you should find around -149.9046

View File

@ -54,14 +54,24 @@ subroutine run
call write_time(6)
call write_int(6,iteration,'CAS-SCF iteration = ')
call write_double(6,energy,'CAS-SCF energy = ')
call write_double(6,energy,'State-average CAS-SCF energy = ')
! if(n_states == 1)then
! call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
! call ezfio_get_casscf_cipsi_energy(PT2)
double precision :: delta_E_istate, e_av
e_av = 0.d0
do istate=1,N_states
call write_double(6,E_PT2(istate),'E + PT2 energy = ')
call write_double(6,PT2(istate),' PT2 = ')
e_av += state_average_weight(istate) * Ev(istate)
if(istate.gt.1)then
delta_E_istate = E_PT2(istate) - E_PT2(1)
write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2 = ',delta_E_istate
endif
write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate)
write(*,'(A6,I2,A18,F16.10)')'state ',istate,' PT2 energy = ',PT2(istate)
! call write_double(6,E_PT2(istate),'E + PT2 energy = ')
! call write_double(6,PT2(istate),' PT2 = ')
enddo
call write_double(6,e_av,'State-average CAS-SCF energy bis = ')
call write_double(6,pt2_max,' PT2_MAX = ')
! endif
@ -99,8 +109,8 @@ subroutine run
mo_coef = NewOrbs
mo_occ = occnum
call save_mos
if(.not.converged)then
call save_mos
iteration += 1
if(norm_grad_vec2.gt.0.01d0)then
N_det = N_states

View File

@ -0,0 +1,66 @@
# This is an example for MULTI STATE CALCULATION STATE AVERAGE CASSCF
# We will compute 3 states on the O2 molecule
# The Ground state and 2 degenerate excited states
# Please follow carefully the tuto :)
##### PREPARING THE EZFIO
# Set the path to your QP2 directory
QP_ROOT=my_fancy_path
source ${QP_ROOT}/quantum_package.rc
# Create the EZFIO folder
qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz_multi_state
# Start with ROHF orbitals
qp run scf # ROHF energy : -149.619992871398
# Freeze the 1s orbitals of the two oxygen
qp set_frozen_core
##### PREPARING THE ORBITALS WITH NATURAL ORBITALS OF A CIS
# Tell that you want 3 states in your WF
qp set determinants n_states 3
# Run a CIS wave function to start your calculation
qp run cis | tee ${EZFIO_FILE}.cis_3_states.out # -149.6652601409258 -149.4714726176746 -149.4686165431939
# Save the STATE AVERAGE natural orbitals for having a balanced description
# This will also order the orbitals according to their occupation number
# Which makes the active space selection easyer !
qp run save_natorb | tee ${EZFIO_FILE}.natorb_3states.out
##### PREPARING A CIS GUESS WITHIN THE ACTIVE SPACE
# Set an active space which has the most of important excitations
# and that maintains symmetry : the ACTIVE ORBITALS are from """6 to 13"""
# YOU FIRST FREEZE THE VIRTUALS THAT ARE NOT IN THE ACTIVE SPACE
# !!!!! WE SET TO "-D" for DELETED !!!!
qp set_mo_class -c "[1-5]" -a "[6-13]" -d "[14-46]"
# You create a guess of CIS type WITHIN THE ACTIVE SPACE
qp run cis | tee ${EZFIO_FILE}.cis_3_states_active_space.out # -149.6515472533511 -149.4622878024821 -149.4622878024817
# You tell to read the WFT stored (i.e. the guess we just created)
qp set determinants read_wf True
##### DOING THE CASSCF
### SETTING PROPERLY THE ACTIVE SPACE FOR CASSCF
# You set the active space WITH THE VIRTUAL ORBITALS !!!
# !!!!! NOW WE SET TO "-v" for VIRTUALS !!!!!
qp set_mo_class -c "[1-5]" -a "[6-13]" -v "[14-46]"
# You tell that it is a small actice space so the CIPSI can take all Slater determinants
qp set casscf_cipsi small_active_space True
# You specify the output file
output=${EZFIO_FILE}.casscf_3states.out
# You run the CASSCF calculation
qp run casscf | tee ${output} # -149.7175867510 -149.5059010227 -149.5059010226
# Some grep in order to get some numbers useful to check convergence
# State average energy
grep "State-average CAS-SCF energy =" $output | cut -d "=" -f 2 > data_e_average
# Delta E anticipated for State-average energy, only usefull to check convergence
grep "Predicted energy improvement =" $output | cut -d "=" -f 2 > data_improve
# Ground state energy
grep "state 1 E + PT2 energy" $output | cut -d "=" -f 2 > data_1
# First excited state energy
grep "state 2 E + PT2 energy" $output | cut -d "=" -f 2 > data_2
# First excitation energy
grep "state 2 Delta E+PT2" $output | cut -d "=" -f 2 > data_delta_E2
# Second excited state energy
grep "state 3 E + PT2 energy" $output | cut -d "=" -f 2 > data_3
# Second excitation energy
grep "state 3 Delta E+PT2" $output | cut -d "=" -f 2 > data_delta_E3

View File

@ -226,27 +226,28 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
end do
! Form the exponential
call exp_matrix_taylor(Tmat,mo_num,Umat,converged)
Tpotmat(:,:)=0.D0
Umat(:,:) =0.D0
do i=1,mo_num
Tpotmat(i,i)=1.D0
Umat(i,i) =1.d0
end do
iter=0
converged=.false.
do while (.not.converged)
iter+=1
f = 1.d0 / dble(iter)
Tpotmat2(:,:) = Tpotmat(:,:) * f
call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, &
Tpotmat2, size(Tpotmat2,1), &
Tmat, size(Tmat,1), 0.d0, &
Tpotmat, size(Tpotmat,1))
Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
end do
! Tpotmat(:,:)=0.D0
! Umat(:,:) =0.D0
! do i=1,mo_num
! Tpotmat(i,i)=1.D0
! Umat(i,i) =1.d0
! end do
! iter=0
! converged=.false.
! do while (.not.converged)
! iter+=1
! f = 1.d0 / dble(iter)
! Tpotmat2(:,:) = Tpotmat(:,:) * f
! call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, &
! Tpotmat2, size(Tpotmat2,1), &
! Tmat, size(Tmat,1), 0.d0, &
! Tpotmat, size(Tpotmat,1))
! Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
!
! converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
! end do
END_PROVIDER

View File

@ -492,3 +492,25 @@ subroutine u_0_H_u_0_two_e(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
deallocate (s_0, v_0)
end
BEGIN_PROVIDER [double precision, psi_energy_two_e_trans, (N_states, N_states)]
implicit none
BEGIN_DOC
! psi_energy_two_e_trans(istate,jstate) = <Psi_istate|W_ee |Psi_jstate>
END_dOC
integer :: i,j,istate,jstate
double precision :: hij, coef_i, coef_j
psi_energy_two_e_trans = 0.d0
do i = 1, N_det
do j = 1, N_det
call i_H_j_two_e(psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
do istate = 1, N_states
coef_i = psi_coef(i,istate)
do jstate = 1, N_states
coef_j = psi_coef(j,jstate)
psi_energy_two_e_trans(jstate,istate) += coef_i * coef_j * hij
enddo
enddo
enddo
enddo
END_PROVIDER

View File

@ -21,3 +21,10 @@ type: logical
doc: If true and N_states > 1, the oscillator strength will be computed
interface: ezfio,provider,ocaml
default: false
[calc_energy_components]
type: logical
doc: If true, the components of the energy (1e, 2e, kinetic) will be computed
interface: ezfio,provider,ocaml
default: false

View File

@ -6,6 +6,11 @@ subroutine print_mol_properties()
! Run the propertie calculations
END_DOC
! Energy components
if (calc_energy_components) then
call print_energy_components
endif
! Electric dipole moment
if (calc_dipole_moment) then
call print_dipole_moment
@ -18,7 +23,7 @@ subroutine print_mol_properties()
! Oscillator strength
if (calc_osc_str .and. N_states > 1) then
call print_oscillator_strength
call print_oscillator_strength
endif
end

View File

@ -0,0 +1,39 @@
BEGIN_PROVIDER [double precision, act_2_rdm_trans_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states,N_states)]
implicit none
BEGIN_DOC
! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) = STATE SPECIFIC physicist notation for 2rdm_trans
!
! \sum_{\sigma,\sigma'}<Psi_{istate}| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi_{jstate}>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
integer :: ispin
double precision :: wall_1, wall_2
! condition for beta/beta spin
print*,''
print*,'Providing act_2_rdm_trans_spin_trace_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_trans_spin_trace_mo'
ispin = 4
act_2_rdm_trans_spin_trace_mo = 0.d0
call wall_time(wall_1)
! if(read_two_body_rdm_trans_spin_trace)then
! print*,'Reading act_2_rdm_trans_spin_trace_mo from disk ...'
! call read_array_two_rdm_trans(n_act_orb,N_states,act_2_rdm_trans_spin_trace_mo,name_file)
! else
call orb_range_2_trans_rdm_openmp(act_2_rdm_trans_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
! endif
! if(write_two_body_rdm_trans_spin_trace)then
! print*,'Writing act_2_rdm_trans_spin_trace_mo on disk ...'
! call write_array_two_rdm_trans(n_act_orb,n_states,act_2_rdm_trans_spin_trace_mo,name_file)
! call ezfio_set_two_body_rdm_trans_io_two_body_rdm_trans_spin_trace("Read")
! endif
act_2_rdm_trans_spin_trace_mo *= 2.d0
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_trans_spin_trace_mo',wall_2 - wall_1
END_PROVIDER

View File

@ -365,3 +365,91 @@ subroutine routine_full_mos
end
subroutine routine_active_only_trans
implicit none
integer :: i,j,k,l,iorb,jorb,korb,lorb,istate,jstate
BEGIN_DOC
! This routine computes the two electron repulsion within the active space using various providers
!
END_DOC
double precision :: vijkl,get_two_e_integral
double precision :: wee_tot(N_states,N_states),rdm_transtot
double precision :: spin_trace
double precision :: accu_tot
wee_tot = 0.d0
iorb = 1
jorb = 1
korb = 1
lorb = 1
vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map)
provide act_2_rdm_trans_spin_trace_mo
i = 1
j = 2
print*,'**************************'
print*,'**************************'
do jstate = 1, N_states
do istate = 1, N_states
!! PURE ACTIVE PART
!!
accu_tot = 0.d0
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_act_orb
jorb = list_act(j)
do k = 1, n_act_orb
korb = list_act(k)
do l = 1, n_act_orb
lorb = list_act(l)
! 1 2 1 2 2 1 2 1
! if(dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate)).gt.1.d-10)then
! print*,'Error in act_2_rdm_trans_spin_trace_mo'
! print*,"dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l) - act_2_rdm_trans_spin_trace_mo(j,i,l,k)).gt.1.d-10"
! print*,i,j,k,l
! print*,act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate),act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate),dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate))
! endif
! 1 2 1 2 1 2 1 2
! if(dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)).gt.1.d-10)then
! print*,'Error in act_2_rdm_trans_spin_trace_mo'
! print*,"dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)).gt.1.d-10"
! print*,i,j,k,l
! print*,act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate),act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate),dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate))
! endif
vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map)
rdm_transtot = act_2_rdm_trans_spin_trace_mo(l,k,j,i,istate,jstate)
wee_tot(istate,jstate) += 0.5d0 * vijkl * rdm_transtot
enddo
enddo
enddo
enddo
print*,''
print*,''
print*,'Active space only energy for state ',istate,jstate
print*,'wee_tot = ',wee_tot(istate,jstate)
print*,'Full energy '
print*,'psi_energy_two_e(istate,jstate)= ',psi_energy_two_e_trans(istate,jstate)
print*,'--------------------------'
enddo
enddo
print*,'Wee from DM '
do istate = 1,N_states
write(*,'(100(F16.10,X))')wee_tot(1:N_states,istate)
enddo
print*,'Wee from Psi det'
do istate = 1,N_states
write(*,'(100(F16.10,X))')psi_energy_two_e_trans(1:N_states,istate)
enddo
end

View File

@ -31,3 +31,37 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file)
close(unit=i_unit_output)
end
subroutine write_array_two_trans_rdm(n_orb,nstates,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb,nstates
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates,nstates)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
call lock_io()
write(i_unit_output)array_tmp
call unlock_io()
close(unit=i_unit_output)
end
subroutine read_array_two_trans_rdm(n_orb,nstates,array_tmp,name_file)
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb,nstates
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states,nstates)
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
call lock_io()
read(i_unit_output)array_tmp
call unlock_io()
close(unit=i_unit_output)
end

View File

@ -4,5 +4,6 @@ program test_2_rdm
touch read_wf
call routine_active_only
call routine_full_mos
call routine_active_only_trans
end

View File

@ -0,0 +1,585 @@
subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! if ispin == 1 :: alpha/alpha 2_rdm
! == 2 :: beta /beta 2_rdm
! == 3 :: alpha/beta + beta/alpha 2trans_rdm
! == 4 :: spin traced 2_rdm :: aa + bb + ab + ba
!
! notice that here it is the TRANSITION RDM THAT IS COMPUTED
!
! THE DIAGONAL PART IS THE USUAL ONE FOR A GIVEN STATE
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st)
double precision, intent(in) :: u_0(sze,N_st)
integer :: k
double precision, allocatable :: u_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
PROVIDE mo_two_e_integrals_in_map
allocate(u_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1)
deallocate(u_t)
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes two-trans_rdm
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st)
double precision, intent(in) :: u_t(N_st,N_det)
integer :: k
PROVIDE N_int
select case (N_int)
case (1)
call orb_range_2_trans_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case (2)
call orb_range_2_trans_rdm_openmp_work_2(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case (3)
call orb_range_2_trans_rdm_openmp_work_3(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case (4)
call orb_range_2_trans_rdm_openmp_work_4(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case default
call orb_range_2_trans_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
use omp_lib
implicit none
BEGIN_DOC
! Computes the two trans_rdm for the N_st vectors |u_t>
! if ispin == 1 :: alpha/alpha 2trans_rdm
! == 2 :: beta /beta 2trans_rdm
! == 3 :: alpha/beta 2trans_rdm
! == 4 :: spin traced 2trans_rdm :: aa + bb + 0.5 (ab + ba))
! The 2trans_rdm will be computed only on the list of orbitals list_orb, which contains norb
! Default should be 1,N_det,0,1 for istart,iend,ishift,istep
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st)
integer(omp_lock_kind) :: lock_2trans_rdm
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b
integer :: krow, kcol
integer :: lrow, lcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
integer(bit_kind) :: orb_bitmask($N_int)
integer :: list_orb_reverse(mo_num)
integer, allocatable :: keys(:,:)
double precision, allocatable :: values(:,:,:)
integer :: nkeys,sze_buff
integer :: ll
alpha_alpha = .False.
beta_beta = .False.
alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
else
print*,'Wrong parameter for ispin in general_2_trans_rdm_state_av_openmp_work'
print*,'ispin = ',ispin
stop
endif
PROVIDE N_int
call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60
list_orb_reverse = -1000
do i = 1, norb
list_orb_reverse(list_orb(i)) = i
enddo
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
call omp_init_lock(lock_2trans_rdm)
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson elec_alpha_num
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2trans_rdm,&
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int,&
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_st, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, &
!$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, &
!$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, &
!$OMP lcol, lrow, l_a, l_b, &
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP n_singles_b, nkeys, keys, values)
! Alpha/Beta double excitations
! =============================
nkeys = 0
allocate( keys(4,sze_buff), values(n_st,n_st,sze_buff))
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
if(alpha_beta.or.spin_trace)then
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
! print*,'nkeys before = ',nkeys
do ll = 1, N_states
do l= 1, N_states
c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a)
enddo
enddo
if(alpha_beta)then
! only ONE contribution
if (nkeys+1 .ge. sze_buff) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
else if (spin_trace)then
! TWO contributions
if (nkeys+2 .ge. sze_buff) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
endif
call orb_range_off_diag_double_to_all_states_ab_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
enddo
endif
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha exitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
do ll= 1, N_states
do l= 1, N_states
c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a)
enddo
enddo
if(alpha_beta.or.spin_trace.or.alpha_alpha)then
! increment the alpha/beta part for single excitations
if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! increment the alpha/alpha part for single excitations
if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
endif
enddo
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
! Compute Hij for all alpha doubles
! ----------------------------------
if(alpha_alpha.or.spin_trace)then
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
do ll= 1, N_states
do l= 1, N_states
c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a)
enddo
enddo
if (nkeys+4 .ge. sze_buff) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
call orb_range_off_diag_double_to_all_states_aa_trans_rdm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
enddo
endif
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
l_a = psi_bilinear_matrix_transp_order(l_b)
do ll= 1, N_states
do l= 1, N_states
c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a)
enddo
enddo
if(alpha_beta.or.spin_trace.or.beta_beta)then
! increment the alpha/beta part for single excitations
if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! increment the beta /beta part for single excitations
if (nkeys+4 * elec_alpha_num .ge. sze_buff) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_bb_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
endif
enddo
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
! Compute Hij for all beta doubles
! ----------------------------------
if(beta_beta.or.spin_trace)then
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
l_a = psi_bilinear_matrix_transp_order(l_b)
do ll= 1, N_states
do l= 1, N_states
c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a)
enddo
enddo
if (nkeys+4 .ge. sze_buff) then
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
endif
call orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! print*,'to do orb_range_off_diag_double_to_2_trans_rdm_bb_dm_buffer'
ASSERT (l_a <= N_det)
enddo
endif
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
double precision, external :: diag_wee_mat_elem, diag_S_mat_elem
double precision :: c_1(N_states,N_states)
do ll = 1, N_states
do l = 1, N_states
c_1(l,ll) = u_t(ll,k_a) * u_t(l,k_a)
enddo
enddo
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
call orb_range_diag_to_all_states_2_rdm_trans_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm)
nkeys = 0
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE
subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
use omp_lib
implicit none
integer, intent(in) :: n_st,nkeys,dim1
integer, intent(in) :: keys(4,nkeys)
double precision, intent(in) :: values(n_st,n_st,nkeys)
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st,n_st)
integer(omp_lock_kind),intent(inout):: lock_2rdm
integer :: i,h1,h2,p1,p2,istate,jstate
call omp_set_lock(lock_2rdm)
! print*,'*************'
! print*,'updating'
! print*,'nkeys',nkeys
do i = 1, nkeys
h1 = keys(1,i)
h2 = keys(2,i)
p1 = keys(3,i)
p2 = keys(4,i)
do jstate = 1, N_st
do istate = 1, N_st
!! print*,h1,h2,p1,p2,values(istate,i)
big_array(h1,h2,p1,p2,istate,jstate) += values(istate,jstate,i)
enddo
enddo
enddo
call omp_unset_lock(lock_2rdm)
end

File diff suppressed because it is too large Load Diff