From 5aab702257bf7213616570a2b177e1c0c766b6d9 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 21 Feb 2023 18:20:15 +0100 Subject: [PATCH 01/16] removed stupid print in dft_utils_in_r/dm_in_r.irp.f --- external/qp2-dependencies | 2 +- src/dft_utils_in_r/dm_in_r.irp.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/external/qp2-dependencies b/external/qp2-dependencies index b8cd5815..f40bde09 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit b8cd5815bce14c9b880e3c5ea3d5fc2652f5955c +Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8 diff --git a/src/dft_utils_in_r/dm_in_r.irp.f b/src/dft_utils_in_r/dm_in_r.irp.f index feb174fd..e4e56ebf 100644 --- a/src/dft_utils_in_r/dm_in_r.irp.f +++ b/src/dft_utils_in_r/dm_in_r.irp.f @@ -82,7 +82,7 @@ enddo enddo !$OMP END PARALLEL DO - print*,'density and gradients provided' +! print*,'density and gradients provided' END_PROVIDER From 0d6de7ce3d78101b50bf2729ce63e3661751691b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 22 Feb 2023 09:51:08 +0100 Subject: [PATCH 02/16] Created codemeta.json --- codemeta.json | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 codemeta.json diff --git a/codemeta.json b/codemeta.json new file mode 100644 index 00000000..9e63b5bc --- /dev/null +++ b/codemeta.json @@ -0,0 +1,66 @@ +{ + "@context": "https://doi.org/10.5063/schema/codemeta-2.0", + "@type": "SoftwareSourceCode", + "license": "https://spdx.org/licenses/AGPL-3.0", + "codeRepository": "https://github.com/QuantumPackage/qp2", + "dateCreated": "2014-04-01", + "datePublished": "2019-06-11", + "dateModified": "2020-02-21", + "downloadUrl": "https://github.com/QuantumPackage/qp2/releases/tag/2.1.2", + "issueTracker": "https://github.com/QuantumPackage/qp2/issues", + "name": "Quantum Package", + "version": "2.1.2", + "identifier": "https://doi.org/10.5281/zenodo.3677565", + "description": "Programming environment for wave function methods", + "applicationCategory": "Quantum Chemistry", + "funding": "ERC_863481, CoE_952165", + "developmentStatus": "active", + "referencePublication": "https://doi.org/10.1021/acs.jctc.9b00176", + "funder": { + "@type": "Organization", + "name": "CNRS" + }, + "keywords": [ + "selected configuration interaction", + "CIPSI" + ], + "programmingLanguage": [ + "Fortran", + "IRPF90", + "OCaml", + "Python", + "C" + ], + "operatingSystem": [ + "Linux" + ], + "softwareRequirements": [ + "ZeroMQ" + ], + "author": [ + { + "@type": "Person", + "@id": "https://orcid.org/0000-0003-4955-7136", + "givenName": "Scemama", + "familyName": "Anthony", + "email": "scemama@irsamc.ups-tlse.fr", + "affiliation": { + "@type": "Organization", + "name": "Laboratoire de chimie et physique quantiques, Toulouse, CNRS" + } + } + ], + "contributor": [ + { + "@type": "Person", + "givenName": "Emmanuel", + "familyName": "Giner", + "email": "eginer@lct.jussieu.fr", + "affiliation": { + "@type": "Organization", + "name": "Laboratoire de Chimie Theorique, Paris, CNRS" + } + } + ] +} + From 656709b7c12de4d5ea7c47d06621fb1a20f6cd4d Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 23 Feb 2023 16:12:00 +0100 Subject: [PATCH 03/16] added tc spin density --- src/tc_bi_ortho/spin_mulliken.irp.f | 106 ++++++++++++++++++++ src/tc_bi_ortho/tc_natorb.irp.f | 2 +- src/tc_bi_ortho/tc_prop.irp.f | 145 +++++++++++++++++++++------- 3 files changed, 218 insertions(+), 35 deletions(-) create mode 100644 src/tc_bi_ortho/spin_mulliken.irp.f diff --git a/src/tc_bi_ortho/spin_mulliken.irp.f b/src/tc_bi_ortho/spin_mulliken.irp.f new file mode 100644 index 00000000..922cc1b9 --- /dev/null +++ b/src/tc_bi_ortho/spin_mulliken.irp.f @@ -0,0 +1,106 @@ + +BEGIN_PROVIDER [double precision, tc_spin_population, (ao_num,ao_num,N_states)] + implicit none + integer :: i,j,istate + BEGIN_DOC +! spin population on the ao basis : +! tc_spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * + END_DOC + tc_spin_population = 0.d0 + do istate = 1, N_states + do i = 1, ao_num + do j = 1, ao_num + tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i) + enddo + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum, (0:ao_l_max,N_states)] +&BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum_per_atom, (0:ao_l_max,nucl_num,N_states)] + implicit none + integer :: i,istate + double precision :: accu + tc_spin_population_angular_momentum = 0.d0 + tc_spin_population_angular_momentum_per_atom = 0.d0 + do istate = 1, N_states + do i = 1, ao_num + tc_spin_population_angular_momentum(ao_l(i),istate) += tc_spin_gross_orbital_product(i,istate) + tc_spin_population_angular_momentum_per_atom(ao_l(i),ao_nucl(i),istate) += tc_spin_gross_orbital_product(i,istate) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, tc_spin_gross_orbital_product, (ao_num,N_states)] + implicit none + tc_spin_gross_orbital_product = 0.d0 + integer :: i,j,istate + BEGIN_DOC +! gross orbital product for the spin population + END_DOC + do istate = 1, N_states + do i = 1, ao_num + do j = 1, ao_num + tc_spin_gross_orbital_product(i,istate) += tc_spin_population(j,i,istate) + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [double precision, tc_mulliken_spin_densities, (nucl_num,N_states)] + implicit none + integer :: i,j,istate + BEGIN_DOC +!ATOMIC SPIN POPULATION (ALPHA MINUS BETA) + END_DOC + tc_mulliken_spin_densities = 0.d0 + do istate = 1, N_states + do i = 1, ao_num + tc_mulliken_spin_densities(ao_nucl(i),istate) += tc_spin_gross_orbital_product(i,istate) + enddo + enddo + +END_PROVIDER + +subroutine tc_print_mulliken_sd + implicit none + double precision :: accu + integer :: i + integer :: j + print*,'Mulliken spin densities' + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),tc_mulliken_spin_densities(i,1) + accu += tc_mulliken_spin_densities(i,1) + enddo + print*,'Sum of Mulliken SD = ',accu + print*,'AO SPIN POPULATIONS' + accu = 0.d0 + do i = 1, ao_num + accu += tc_spin_gross_orbital_product(i,1) + write(*,'(1X,I3,1X,A4,1X,I2,1X,A4,1X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_character(ao_l(i))),tc_spin_gross_orbital_product(i,1) + enddo + print*,'sum = ',accu + accu = 0.d0 + print*,'Angular momentum analysis' + do i = 0, ao_l_max + accu += tc_spin_population_angular_momentum(i,1) + print*,' ',trim(l_to_character(i)),tc_spin_population_angular_momentum(i,1) + print*,'sum = ',accu + enddo + print*,'Angular momentum analysis per atom' + print*,'Angular momentum analysis' + do j = 1,nucl_num + accu = 0.d0 + do i = 0, ao_l_max + accu += tc_spin_population_angular_momentum_per_atom(i,j,1) + write(*,'(1X,I3,1X,A4,1X,A4,1X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_character(i)),tc_spin_population_angular_momentum_per_atom(i,j,1) + print*,'sum = ',accu + enddo + enddo + +end + diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index 33410570..b7e5ae81 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -21,7 +21,7 @@ allocate(dm_tmp(mo_num,mo_num), fock_diag(mo_num)) - dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1) + dm_tmp(1:mo_num,1:mo_num) = -tc_transition_matrix_mo(1:mo_num,1:mo_num,1,1) print *, ' dm_tmp' do i = 1, mo_num diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f index c7f6c986..5bb0e2c0 100644 --- a/src/tc_bi_ortho/tc_prop.irp.f +++ b/src/tc_bi_ortho/tc_prop.irp.f @@ -1,8 +1,11 @@ -BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_states,N_states) ] + BEGIN_PROVIDER [ double precision, tc_transition_matrix_mo_beta, (mo_num, mo_num,N_states,N_states) ] +&BEGIN_PROVIDER [ double precision, tc_transition_matrix_mo_alpha, (mo_num, mo_num,N_states,N_states) ] implicit none BEGIN_DOC - ! tc_transition_matrix(p,h,istate,jstate) = + ! tc_transition_matrix_mo_alpha(p,h,istate,jstate) = + ! + ! tc_transition_matrix_mo_beta(p,h,istate,jstate) = ! ! where are the left/right eigenvectors on a bi-ortho basis END_DOC @@ -11,43 +14,65 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state integer, allocatable :: occ(:,:) integer :: n_occ_ab(2),degree,exc(0:2,2,2) allocate(occ(N_int*bit_kind_size,2)) - tc_transition_matrix = 0.d0 - do istate = 1, N_states - do jstate = 1, N_states + tc_transition_matrix_mo_alpha = 0.d0 + tc_transition_matrix_mo_beta = 0.d0 do i = 1, N_det do j = 1, N_det call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) - if(degree.gt.1)then - cycle - else if (degree == 0)then - call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) - do p = 1, n_occ_ab(1) ! browsing the alpha electrons - m = occ(p,1) - tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) - enddo - do p = 1, n_occ_ab(2) ! browsing the beta electrons - m = occ(p,1) - tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) - enddo - else - call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int) - if (exc(0,1,1) == 1) then - ! Single alpha - h = exc(1,1,1) ! hole in psi_det(1,1,j) - p = exc(1,2,1) ! particle in psi_det(1,1,j) - else - ! Single beta - h = exc(1,1,2) ! hole in psi_det(1,1,j) - p = exc(1,2,2) ! particle in psi_det(1,1,j) - endif - tc_transition_matrix(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) - endif + if(degree.gt.1)cycle + do istate = 1, N_states + do jstate = 1, N_states + if (degree == 0)then + call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + do p = 1, n_occ_ab(1) ! browsing the alpha electrons + m = occ(p,1) + tc_transition_matrix_mo_alpha(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + enddo + do p = 1, n_occ_ab(2) ! browsing the beta electrons + m = occ(p,1) + tc_transition_matrix_mo_beta(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + enddo + else + call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Single alpha + h = exc(1,1,1) ! hole in psi_det(1,1,j) + p = exc(1,2,1) ! particle in psi_det(1,1,j) + tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + else + ! Single beta + h = exc(1,1,2) ! hole in psi_det(1,1,j) + p = exc(1,2,2) ! particle in psi_det(1,1,j) + tc_transition_matrix_mo_beta(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + endif + endif + enddo enddo - enddo enddo enddo END_PROVIDER + BEGIN_PROVIDER [double precision, tc_transition_matrix_mo, (mo_num, mo_num,N_states,N_states) ] + implicit none + BEGIN_DOC + ! tc_transition_matrix_mo(p,h,istate,jstate) = \sum_{sigma=alpha,beta} + ! + ! where are the left/right eigenvectors on a bi-ortho basis + END_DOC + tc_transition_matrix_mo = tc_transition_matrix_mo_beta + tc_transition_matrix_mo_alpha + END_PROVIDER + + + BEGIN_PROVIDER [double precision, tc_spin_transition_matrix_mo, (mo_num, mo_num,N_states,N_states) ] + implicit none + BEGIN_DOC + ! tc_spin_transition_matrix_mo = tc_transition_matrix_mo_alpha - tc_transition_matrix_mo_beta + ! + ! where are the left/right eigenvectors on a bi-ortho basis + END_DOC + tc_spin_transition_matrix_mo = tc_transition_matrix_mo_alpha - tc_transition_matrix_mo_beta + END_PROVIDER + BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)] implicit none @@ -57,9 +82,9 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state do istate = 1, N_states do i = 1, mo_num do j = 1, mo_num - tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i) - tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i) - tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i) + tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i) + tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i) + tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i) enddo enddo enddo @@ -78,3 +103,55 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state enddo END_PROVIDER + + BEGIN_PROVIDER [ double precision, tc_transition_matrix_ao, (ao_num, ao_num,N_states,N_states) ] + implicit none + BEGIN_DOC +! tc_transition_matrix(p,h,istate,jstate) in the AO basis + END_DOC + integer :: i,j,k,l + double precision :: dm_mo + tc_transition_matrix_ao = 0.d0 + integer :: istate,jstate + do istate = 1, N_states + do jstate = 1, N_states + do i = 1, mo_num + do j = 1, mo_num + dm_mo = tc_transition_matrix_mo(j,i,jstate,istate) + do k = 1, ao_num + do l = 1, ao_num + tc_transition_matrix_ao(l,k,jstate,istate) += mo_l_coef(l,j) * mo_r_coef(k,i) * dm_mo + enddo + enddo + enddo + enddo + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [ double precision, tc_spin_transition_matrix_ao, (ao_num, ao_num,N_states,N_states) ] + implicit none + BEGIN_DOC +! tc_spin_transition_matrix_ao(p,h,istate,jstate) in the AO basis + END_DOC + integer :: i,j,k,l + double precision :: dm_mo + tc_spin_transition_matrix_ao = 0.d0 + integer :: istate,jstate + do istate = 1, N_states + do jstate = 1, N_states + do i = 1, mo_num + do j = 1, mo_num + dm_mo = tc_spin_transition_matrix_mo(j,i,jstate,istate) + do k = 1, ao_num + do l = 1, ao_num + tc_spin_transition_matrix_ao(l,k,jstate,istate) += mo_l_coef(l,j) * mo_r_coef(k,i) * dm_mo + enddo + enddo + enddo + enddo + enddo + enddo + + END_PROVIDER From ef3a52f54d5ef00d4abbc310b7e5beec1851b7e2 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 23 Feb 2023 16:38:40 +0100 Subject: [PATCH 04/16] added write_pt_charges.py --- src/hartree_fock/10.hf.bats | 20 +++++---- src/nuclei/write_pt_charges.py | 80 ++++++++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 9 deletions(-) create mode 100644 src/nuclei/write_pt_charges.py diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index 20b59603..df566032 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -34,23 +34,28 @@ cat > hcn_charges.xyz << EOF 0.5 -2.0 0.0 0.0 EOF -rm -rf hcn.ezfio -qp create_ezfio -b def2-svp hcn.xyz +EZFIO=hcn_pt_charges +rm -rf $EZFIO +qp create_ezfio -b def2-svp hcn.xyz -o $EZFIO qp run scf -mv hcn_charges.xyz hcn.ezfio_point_charges.xyz -python write_pt_charges.py hcn.ezfio +mv hcn_charges.xyz ${EZFIO}_point_charges.xyz +python write_pt_charges.py ${EZFIO} qp set nuclei point_charges True -qp run scf | tee hcn.ezfio.pt_charges.out +qp run scf | tee ${EZFIO}.pt_charges.out energy="$(ezfio get hartree_fock energy)" -rm -rf hcn.ezfio good=-92.76613324421798 eq $energy $good $thresh +rm -rf $EZFIO } @test "point charges" { run_pt_charges } +@test "HCN" { # 7.792500 8.51926s + run hcn.ezfio -92.88717500035233 +} + @test "B-B" { # 3s run b2_stretched.ezfio -48.9950585434279 } @@ -124,9 +129,6 @@ good=-92.76613324421798 run ch4.ezfio -40.19961807784367 } -@test "HCN" { # 7.792500 8.51926s - run hcn.ezfio -92.88717500035233 -} @test "N2" { # 8.648100 13.754s run n2.ezfio -108.9834897852979 diff --git a/src/nuclei/write_pt_charges.py b/src/nuclei/write_pt_charges.py new file mode 100644 index 00000000..6dbcd5b8 --- /dev/null +++ b/src/nuclei/write_pt_charges.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python +import os +import sys + +# First argument is the EZFIO file +# It reads a file EZFIO_point_charges.xyz written in this way: +# charge x y z (Angstrom) +# for all charges + + +def zip_in_ezfio(ezfio,tmp): + tmpzip=tmp+".gz" + cmdzip="gzip -c "+tmp+" > "+tmpzip + os.system(cmdzip) + os.system("rm "+tmp) + cmdmv="mv "+tmpzip+" "+EZFIO+"/nuclei/"+tmpzip + os.system(cmdmv) + +def mv_in_ezfio(ezfio,tmp): + cmdmv="mv "+tmp+" "+EZFIO+"/nuclei/"+tmp + os.system(cmdmv) + + +# Getting the EZFIO +EZFIO=sys.argv[1] +EZFIO=EZFIO.replace("/", "") +print(EZFIO) + +# Reading the point charges and convert the Angstrom geometry in Bohr for QP +f = open(EZFIO+'_point_charges.xyz','r') +lines = f.readlines() +convert_angs_to_bohr=1.8897259885789233 + +n_charges=0 +coord_x=[] +coord_y=[] +coord_z=[] +charges=[] +for line in lines: + data = line.split() + if(len(data)>0): + n_charges += 1 + charges.append(str(data[0])) + coord_x.append(str(convert_angs_to_bohr*float(data[1]))) + coord_y.append(str(convert_angs_to_bohr*float(data[2]))) + coord_z.append(str(convert_angs_to_bohr*float(data[3]))) + +# Write the file containing the number of charges and set in EZFIO folder +tmp="n_pts_charge" +fncharges = open(tmp,'w') +fncharges.write(" "+str(n_charges)+'\n') +fncharges.close() +mv_in_ezfio(EZFIO,tmp) + +# Write the file containing the charges and set in EZFIO folder +tmp="pts_charge_z" +fcharges = open(tmp,'w') +fcharges.write(" 1\n") +fcharges.write(" "+str(n_charges)+'\n') +for i in range(n_charges): + fcharges.write(charges[i]+'\n') +fcharges.close() +zip_in_ezfio(EZFIO,tmp) + +# Write the file containing the charge coordinates and set in EZFIO folder +tmp="pts_charge_coord" +fcoord = open(tmp,'w') +fcoord.write(" 2\n") +fcoord.write(" "+str(n_charges)+' 3\n') +#fcoord.write(" "+' 3 '+str(n_charges)+' \n') +for i in range(n_charges): + fcoord.write(' '+coord_x[i]+'\n') +for i in range(n_charges): + fcoord.write(' '+coord_y[i]+'\n') +for i in range(n_charges): + fcoord.write(' '+coord_z[i]+'\n') +fcoord.close() +zip_in_ezfio(EZFIO,tmp) + + From c4445d9a61911d30b5870e7c259f3cee85ccaccd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 24 Feb 2023 16:33:28 +0100 Subject: [PATCH 05/16] Fix qp2_dependencies branch --- configure | 2 ++ external/qp2-dependencies | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/configure b/configure index ebfb70b4..5b50d0d7 100755 --- a/configure +++ b/configure @@ -19,6 +19,8 @@ git submodule update # Update ARM or x86 dependencies ARCHITECTURE=$(uname -m) cd ${QP_ROOT}/external/qp2-dependencies +git checkout master +git pull echo "Architecture: $ARCHITECTURE" cd ${QP_ROOT} diff --git a/external/qp2-dependencies b/external/qp2-dependencies index f40bde09..b8cd5815 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8 +Subproject commit b8cd5815bce14c9b880e3c5ea3d5fc2652f5955c From 5df44b6b562e03fb38c2b234230438220daf6867 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 24 Feb 2023 18:05:45 +0100 Subject: [PATCH 06/16] Removed Cryptokit --- INSTALL.rst | 2 +- external/qp2-dependencies | 2 +- ocaml/Makefile | 2 +- ocaml/To_md5.ml | 4 ++-- ocaml/_tags | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/INSTALL.rst b/INSTALL.rst index e37d31eb..e9f4aedb 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -316,7 +316,7 @@ OCaml .. code:: bash - opam install ocamlbuild cryptokit zmq sexplib ppx_sexp_conv ppx_deriving getopt + opam install ocamlbuild zmq sexplib ppx_sexp_conv ppx_deriving getopt Docopt diff --git a/external/qp2-dependencies b/external/qp2-dependencies index b8cd5815..ce14f57b 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit b8cd5815bce14c9b880e3c5ea3d5fc2652f5955c +Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4 diff --git a/ocaml/Makefile b/ocaml/Makefile index 8853a991..c03be131 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -43,7 +43,7 @@ $(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.m $(QP_ROOT)/ocaml/element_create_db.byte external_libs: - opam install cryptokit sexplib + opam install sexplib qpackage.odocl: $(MLIFILES) ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl diff --git a/ocaml/To_md5.ml b/ocaml/To_md5.ml index bc6608f9..1485678c 100644 --- a/ocaml/To_md5.ml +++ b/ocaml/To_md5.ml @@ -4,8 +4,8 @@ open Sexplib let to_md5 sexp_of_t t = sexp_of_t t |> Sexp.to_string - |> Cryptokit.hash_string (Cryptokit.Hash.md5 ()) - |> Cryptokit.transform_string (Cryptokit.Hexa.encode ()) + |> Digest.string + |> Digest.to_hex |> MD5.of_string ;; diff --git a/ocaml/_tags b/ocaml/_tags index 55b1c681..0ff23d59 100644 --- a/ocaml/_tags +++ b/ocaml/_tags @@ -1,4 +1,4 @@ -true: package(cryptokit,zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt) +true: package(zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt) true: thread false: profile <*byte> : linkdep(c_bindings.o), custom From 4f071a59fb367f6248aafb4962cceaffbecbe14c Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 24 Feb 2023 21:25:36 +0100 Subject: [PATCH 07/16] added tc_spin_density --- src/tc_bi_ortho/12.tc_bi_ortho.bats | 49 ----------------------------- 1 file changed, 49 deletions(-) delete mode 100644 src/tc_bi_ortho/12.tc_bi_ortho.bats diff --git a/src/tc_bi_ortho/12.tc_bi_ortho.bats b/src/tc_bi_ortho/12.tc_bi_ortho.bats deleted file mode 100644 index f5b9d8c0..00000000 --- a/src/tc_bi_ortho/12.tc_bi_ortho.bats +++ /dev/null @@ -1,49 +0,0 @@ -#!/usr/bin/env bats - -source $QP_ROOT/tests/bats/common.bats.sh -source $QP_ROOT/quantum_package.rc - - -function run_Ne() { - qp set_file Ne_tc_scf - qp run cisd - qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out - eref=-128.77020441279302 - energy="$(grep "eigval_right_tc_bi_orth =" Ne_tc_scf.cisd_tc_bi_ortho.out)" - eq $energy $eref 1e-6 -} - - -@test "Ne" { - run_Ne -} - - -function run_C() { - qp set_file C_tc_scf - qp run cisd - qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out - eref=-37.757536149952514 - energy="$(grep "eigval_right_tc_bi_orth =" C_tc_scf.cisd_tc_bi_ortho.out)" - eq $energy $eref 1e-6 -} - - -@test "C" { - run_C -} - -function run_O() { - qp set_file C_tc_scf - qp run cisd - qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out - eref=-74.908518517716161 - energy="$(grep "eigval_right_tc_bi_orth =" O_tc_scf.cisd_tc_bi_ortho.out)" - eq $energy $eref 1e-6 -} - - -@test "O" { - run_O -} - From 274e903d3c5543a6d5446fb8b52271c7ff896246 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 24 Feb 2023 21:38:09 +0100 Subject: [PATCH 08/16] added spin density --- src/tc_bi_ortho/31.tc_bi_ortho.bats | 49 ++++++++++++++++++++++++ src/tc_bi_ortho/print_tc_spin_dens.irp.f | 16 ++++++++ src/tc_bi_ortho/spin_mulliken.irp.f | 14 +++++-- src/tc_bi_ortho/test_spin_dens.irp.f | 30 +++++++++++++++ src/tc_keywords/EZFIO.cfg | 6 +++ 5 files changed, 112 insertions(+), 3 deletions(-) create mode 100644 src/tc_bi_ortho/31.tc_bi_ortho.bats create mode 100644 src/tc_bi_ortho/print_tc_spin_dens.irp.f create mode 100644 src/tc_bi_ortho/test_spin_dens.irp.f diff --git a/src/tc_bi_ortho/31.tc_bi_ortho.bats b/src/tc_bi_ortho/31.tc_bi_ortho.bats new file mode 100644 index 00000000..f5b9d8c0 --- /dev/null +++ b/src/tc_bi_ortho/31.tc_bi_ortho.bats @@ -0,0 +1,49 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh +source $QP_ROOT/quantum_package.rc + + +function run_Ne() { + qp set_file Ne_tc_scf + qp run cisd + qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out + eref=-128.77020441279302 + energy="$(grep "eigval_right_tc_bi_orth =" Ne_tc_scf.cisd_tc_bi_ortho.out)" + eq $energy $eref 1e-6 +} + + +@test "Ne" { + run_Ne +} + + +function run_C() { + qp set_file C_tc_scf + qp run cisd + qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out + eref=-37.757536149952514 + energy="$(grep "eigval_right_tc_bi_orth =" C_tc_scf.cisd_tc_bi_ortho.out)" + eq $energy $eref 1e-6 +} + + +@test "C" { + run_C +} + +function run_O() { + qp set_file C_tc_scf + qp run cisd + qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out + eref=-74.908518517716161 + energy="$(grep "eigval_right_tc_bi_orth =" O_tc_scf.cisd_tc_bi_ortho.out)" + eq $energy $eref 1e-6 +} + + +@test "O" { + run_O +} + diff --git a/src/tc_bi_ortho/print_tc_spin_dens.irp.f b/src/tc_bi_ortho/print_tc_spin_dens.irp.f new file mode 100644 index 00000000..8308140d --- /dev/null +++ b/src/tc_bi_ortho/print_tc_spin_dens.irp.f @@ -0,0 +1,16 @@ +program test_spin_dens + implicit none + BEGIN_DOC +! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call tc_print_mulliken_sd +! call test + +end diff --git a/src/tc_bi_ortho/spin_mulliken.irp.f b/src/tc_bi_ortho/spin_mulliken.irp.f index 922cc1b9..e225d299 100644 --- a/src/tc_bi_ortho/spin_mulliken.irp.f +++ b/src/tc_bi_ortho/spin_mulliken.irp.f @@ -7,13 +7,21 @@ BEGIN_PROVIDER [double precision, tc_spin_population, (ao_num,ao_num,N_states)] ! tc_spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * END_DOC tc_spin_population = 0.d0 - do istate = 1, N_states + if(only_spin_tc_right)then do i = 1, ao_num do j = 1, ao_num - tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i) + tc_spin_population(j,i,1) = tc_spin_dens_right_only(j,i) * ao_overlap(j,i) enddo enddo - enddo + else + do istate = 1, N_states + do i = 1, ao_num + do j = 1, ao_num + tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i) + enddo + enddo + enddo + endif END_PROVIDER BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum, (0:ao_l_max,N_states)] diff --git a/src/tc_bi_ortho/test_spin_dens.irp.f b/src/tc_bi_ortho/test_spin_dens.irp.f new file mode 100644 index 00000000..2c2f6e7c --- /dev/null +++ b/src/tc_bi_ortho/test_spin_dens.irp.f @@ -0,0 +1,30 @@ +BEGIN_PROVIDER [ double precision, mo_r_coef_normalized, (ao_num,mo_num) ] + implicit none + integer :: i,j,k + double precision :: norm + do i = 1, mo_num + norm = 0.d0 + do j = 1, ao_num + do k = 1, ao_num + norm += mo_r_coef(k,i) * mo_r_coef(j,i) * ao_overlap(k,j) + enddo + enddo + norm = 1.d0/dsqrt(norm) + do j = 1, ao_num + mo_r_coef_normalized(j,i) = mo_r_coef(j,i) * norm + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, tc_spin_dens_right_only, (ao_num, ao_num)] + implicit none + integer :: i,j,k + tc_spin_dens_right_only = 0.d0 + do i = elec_beta_num+1, elec_alpha_num + do j = 1, ao_num + do k = 1, ao_num + tc_spin_dens_right_only(k,j) += mo_r_coef_normalized(k,i) * mo_r_coef_normalized(j,i) + enddo + enddo + enddo +END_PROVIDER diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 5d5477bc..4a1fcb9f 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -183,3 +183,9 @@ type: integer doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body interface: ezfio,provider,ocaml default: -1 + +[only_spin_tc_right] +type: logical +doc: If |true|, only the right part of WF is used to compute spin dens +interface: ezfio,provider,ocaml +default: False From 2b57c0728228adae01ae341b7068356f47d488dc Mon Sep 17 00:00:00 2001 From: ydamour Date: Sat, 25 Feb 2023 01:37:41 +0100 Subject: [PATCH 09/16] more efficient restore symmetry --- src/utils/linear_algebra.irp.f | 145 ++++++++++++++++++++------------- 1 file changed, 88 insertions(+), 57 deletions(-) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 51df33c5..c02560e3 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1574,79 +1574,110 @@ subroutine nullify_small_elements(m,n,A,LDA,thresh) end subroutine restore_symmetry(m,n,A,LDA,thresh) + implicit none + BEGIN_DOC -! Tries to find the matrix elements that are the same, and sets them -! to the average value. -! If restore_symm is False, only nullify small elements + ! Tries to find the matrix elements that are the same, and sets them + ! to the average value. + ! If restore_symm is False, only nullify small elements END_DOC + integer, intent(in) :: m,n,LDA double precision, intent(inout) :: A(LDA,n) double precision, intent(in) :: thresh - integer :: i,j,k,l - logical, allocatable :: done(:,:) - double precision :: f, g, count, thresh2 + + double precision, allocatable :: copy(:), copy_sign(:) + integer, allocatable :: key(:), ii(:), jj(:) + integer :: sze, pi, pf, idx, i,j,k + double precision :: average, val, thresh2 + thresh2 = dsqrt(thresh) - call nullify_small_elements(m,n,A,LDA,thresh) -! if (.not.restore_symm) then -! return -! endif + sze = m * n - ! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)): - ! - copy all values in a 1D array - ! - sort 1D array - ! - average nearby elements - ! - for all elements, find matching value in the sorted 1D array + allocate(copy(sze),copy_sign(sze),key(sze),ii(sze),jj(sze)) - allocate(done(m,n)) - - do j=1,n - do i=1,m - done(i,j) = A(i,j) == 0.d0 + ! Copy to 1D + !$OMP PARALLEL if (m>100) & + !$OMP SHARED(A,m,n,sze,copy_sign,copy,key,ii,jj) & + !$OMP PRIVATE(i,j,k) & + !$OMP DEFAULT(NONE) + !$OMP DO + do j = 1, n + do i = 1, m + k = i+(j-1)*m + copy(k) = A(i,j) + copy_sign(k) = sign(1.d0,copy(k)) + copy(k) = -dabs(copy(k)) + key(k) = k + ii(k) = i + jj(k) = j enddo enddo + !$OMP END DO + !$OMP END PARALLEL - do j=1,n - do i=1,m - if ( done(i,j) ) cycle - done(i,j) = .True. - count = 1.d0 - f = 1.d0/A(i,j) - do l=1,n - do k=1,m - if ( done(k,l) ) cycle - g = f * A(k,l) - if ( dabs(dabs(g) - 1.d0) < thresh2 ) then - count = count + 1.d0 - if (g>0.d0) then - A(i,j) = A(i,j) + A(k,l) - else - A(i,j) = A(i,j) - A(k,l) - end if - endif - enddo - enddo - if (count > 1.d0) then - A(i,j) = A(i,j) / count - do l=1,n - do k=1,m - if ( done(k,l) ) cycle - g = f * A(k,l) - if ( dabs(dabs(g) - 1.d0) < thresh2 ) then - done(k,l) = .True. - if (g>0.d0) then - A(k,l) = A(i,j) - else - A(k,l) = -A(i,j) - end if - endif - enddo - enddo + ! Sort + call dsort(copy,key,sze) + call iset_order(ii,key,sze) + call iset_order(jj,key,sze) + call dset_order(copy_sign,key,sze) + + !TODO + ! Parallelization with OMP + + ! Symmetrize + i = 1 + do while (i < sze) + pi = i + pf = i + + ! Exit if the remaining elements are below thresh + if (-copy(i) <= thresh) exit + + val = 1d0/copy(i) + do while (dabs(val * copy(pf+1) - 1d0) < thresh2) + pf = pf + 1 + ! if pf == sze, copy(pf+1) will not be valid + if (pf == sze) then + exit endif - enddo + ! if pi and pf are different do the average from pi to pf + if (pf - pi > 0) then + average = 0d0 + do j = pi, pf + average = average + copy(j) + enddo + average = average / (pf-pi+1.d0) + do j = pi, pf + copy(j) = average + enddo + ! Update i + i = pf + endif + + ! Update i + i = i + 1 enddo + copy(i:) = 0.d0 + + !$OMP PARALLEL if (sze>10000) & + !$OMP SHARED(m,sze,copy_sign,copy,key,A,ii,jj) & + !$OMP PRIVATE(i,j,k,idx) & + !$OMP DEFAULT(NONE) + ! copy -> A + !$OMP DO + do k = 1, sze + i = ii(k) + j = jj(k) + A(i,j) = sign(copy(k),copy_sign(k)) + enddo + !$OMP END DO + !$OMP END PARALLEL + + deallocate(copy,copy_sign,key,ii,jj) end From fe27080069d6685ff1ff94fa9ce1d78d6cc853c2 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 15:27:39 +0100 Subject: [PATCH 10/16] beginning to introduce a factor 2 in two-rdm --- external/qp2-dependencies | 2 +- src/basis_correction/TODO | 2 ++ src/basis_correction/print_routine.irp.f | 6 +++--- src/two_body_rdm/example.irp.f | 2 +- src/two_body_rdm/test_2_rdm.irp.f | 2 +- src/two_rdm_routines/davidson_like_2rdm.irp.f | 4 ++-- 6 files changed, 10 insertions(+), 8 deletions(-) diff --git a/external/qp2-dependencies b/external/qp2-dependencies index ce14f57b..f40bde09 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4 +Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8 diff --git a/src/basis_correction/TODO b/src/basis_correction/TODO index e28d593a..36c438c8 100644 --- a/src/basis_correction/TODO +++ b/src/basis_correction/TODO @@ -1 +1,3 @@ change all correlation functionals with the pbe_on_top general +factor 2 in two-rdm involved in: + on-top, mu(r), pbe-on-top, sc_basis_corr and so on diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f index 67c5c6c2..c2558d22 100644 --- a/src/basis_correction/print_routine.irp.f +++ b/src/basis_correction/print_routine.irp.f @@ -18,7 +18,7 @@ subroutine print_basis_correction print*, '' print*, 'For more details look at Journal of Chemical Physics 149, 194301 1-15 (2018) ' print*, ' Journal of Physical Chemistry Letters 10, 2931-2937 (2019) ' - print*, ' ???REF SC?' + print*, ' Journal of Chemical Physics 152, 174104 (2020) ' print*, '****************************************' print*, '****************************************' print*, 'mu_of_r_potential = ',mu_of_r_potential @@ -56,14 +56,14 @@ subroutine print_basis_correction print*,'' print*,'********************************************' print*,'********************************************' - print*,'+) PBE-on-top Ecmd functional : (??????? REF-SCF ??????????)' + print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) ' print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization' do istate = 1, N_states write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate) enddo print*,'' print*,'********************************************' - print*,'+) PBE-on-top no spin polarization Ecmd functional : (??????? REF-SCF ??????????)' + print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)' print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION' do istate = 1, N_states write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index de3d97b9..67de9df4 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -120,7 +120,7 @@ subroutine routine_active_only wee_ab(istate) += vijkl * rdmab wee_aa(istate) += vijkl * rdmaa wee_bb(istate) += vijkl * rdmbb - wee_tot(istate) += vijkl * rdmtot + wee_tot(istate) += vijkl * rdmtot enddo enddo diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f index 123261d8..4eb8f9f0 100644 --- a/src/two_body_rdm/test_2_rdm.irp.f +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -2,7 +2,7 @@ program test_2_rdm implicit none read_wf = .True. touch read_wf - call routine_active_only +! call routine_active_only call routine_full_mos end diff --git a/src/two_rdm_routines/davidson_like_2rdm.irp.f b/src/two_rdm_routines/davidson_like_2rdm.irp.f index 2e5aa4d1..ad7a3b21 100644 --- a/src/two_rdm_routines/davidson_like_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_2rdm.irp.f @@ -4,8 +4,8 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz BEGIN_DOC ! if ispin == 1 :: alpha/alpha 2rdm ! == 2 :: beta /beta 2rdm - ! == 3 :: alpha/beta 2rdm - ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) + ! == 3 :: alpha/beta + beta/alpha 2rdm + ! == 4 :: spin traced 2rdm :: aa + bb + ab + ba ! ! Assumes that the determinants are in psi_det ! From 3ba5d3b540424df9cacfe578d24b209ed93e018e Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 15:45:09 +0100 Subject: [PATCH 11/16] factor two introduced in non active only non state average two-rdm, it works with example.irp.f --- src/basis_correction/TODO | 1 + src/two_body_rdm/act_2_rdm.irp.f | 4 ++ src/two_body_rdm/example.irp.f | 47 ++++++++++++------------ src/two_body_rdm/state_av_act_2rdm.irp.f | 9 +++-- src/two_body_rdm/two_e_dm_mo.irp.f | 3 +- 5 files changed, 37 insertions(+), 27 deletions(-) diff --git a/src/basis_correction/TODO b/src/basis_correction/TODO index 36c438c8..64a6ddeb 100644 --- a/src/basis_correction/TODO +++ b/src/basis_correction/TODO @@ -1,3 +1,4 @@ change all correlation functionals with the pbe_on_top general factor 2 in two-rdm involved in: on-top, mu(r), pbe-on-top, sc_basis_corr and so on + casscf : state_av_act_2_rdm_spin_trace_mo diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 41b28aea..61ae4e47 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -44,6 +44,7 @@ endif call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1 + act_2_rdm_ab_mo *= 2.d0 END_PROVIDER @@ -84,6 +85,7 @@ call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1 + act_2_rdm_aa_mo *= 2.d0 END_PROVIDER @@ -124,6 +126,7 @@ call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1 + act_2_rdm_bb_mo *= 2.d0 END_PROVIDER BEGIN_PROVIDER [double precision, act_2_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] @@ -161,6 +164,7 @@ call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read") endif + act_2_rdm_spin_trace_mo *= 2.d0 call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1 END_PROVIDER diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 67de9df4..985f6a5d 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -117,10 +117,10 @@ subroutine routine_active_only endif - wee_ab(istate) += vijkl * rdmab - wee_aa(istate) += vijkl * rdmaa - wee_bb(istate) += vijkl * rdmbb - wee_tot(istate) += vijkl * rdmtot + wee_ab(istate) += 0.5d0 * vijkl * rdmab + wee_aa(istate) += 0.5d0 * vijkl * rdmaa + wee_bb(istate) += 0.5d0 * vijkl * rdmbb + wee_tot(istate) += 0.5d0 * vijkl * rdmtot enddo enddo @@ -144,13 +144,13 @@ subroutine routine_active_only print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) print*,'--------------------------' print*,'accu_aa = ',accu_aa - print*,'N_a (N_a-1)/2 = ', elec_alpha_num*(elec_alpha_num-1)*0.5 + print*,'N_a (N_a-1) = ', elec_alpha_num*(elec_alpha_num-1) print*,'accu_bb = ',accu_bb - print*,'N_b (N_b-1)/2 = ', elec_beta_num*(elec_beta_num-1)*0.5 + print*,'2 N_b (N_b-1) = ', elec_beta_num*(elec_beta_num-1)*2 print*,'accu_ab = ',accu_ab print*,'N_a N_b = ', elec_beta_num*elec_alpha_num print*,'accu_tot = ',accu_tot - print*,'Ne(Ne-1)/2 = ',(elec_num-1)*elec_num * 0.5 + print*,'Ne(Ne-1)/2 = ',(elec_num-1)*elec_num enddo wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 @@ -192,10 +192,10 @@ subroutine routine_active_only print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av) endif - wee_aa_st_av += vijkl * rdm_aa_st_av - wee_bb_st_av += vijkl * rdm_bb_st_av - wee_ab_st_av += vijkl * rdm_ab_st_av - wee_tot_st_av += vijkl * rdm_tot_st_av + wee_aa_st_av += 0.5d0 * vijkl * rdm_aa_st_av + wee_bb_st_av += 0.5d0 * vijkl * rdm_bb_st_av + wee_ab_st_av += 0.5d0 * vijkl * rdm_ab_st_av + wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av enddo enddo enddo @@ -279,10 +279,10 @@ subroutine routine_full_mos rdmbb = full_occ_2_rdm_bb_mo(l,k,j,i,istate) rdmtot = full_occ_2_rdm_spin_trace_mo(l,k,j,i,istate) - wee_ab(istate) += vijkl * rdmab - wee_aa(istate) += vijkl * rdmaa - wee_bb(istate) += vijkl * rdmbb - wee_tot(istate)+= vijkl * rdmtot + wee_ab(istate) += 0.5d0 * vijkl * rdmab + wee_aa(istate) += 0.5d0 * vijkl * rdmaa + wee_bb(istate) += 0.5d0 * vijkl * rdmbb + wee_tot(istate)+= 0.5d0 * vijkl * rdmtot enddo enddo aa_norm(istate) += full_occ_2_rdm_aa_mo(j,i,j,i,istate) @@ -310,18 +310,19 @@ subroutine routine_full_mos print*,'Normalization of two-rdms ' print*,'' print*,'aa_norm(istate) = ',aa_norm(istate) - print*,'N_alpha(N_alpha-1)/2 = ',elec_num_tab(1) * (elec_num_tab(1) - 1)/2 + print*,'N_alpha(N_alpha-1) = ',elec_num_tab(1) * (elec_num_tab(1) - 1) print*,'' print*,'bb_norm(istate) = ',bb_norm(istate) - print*,'N_alpha(N_alpha-1)/2 = ',elec_num_tab(2) * (elec_num_tab(2) - 1)/2 + print*,'N_alpha(N_alpha-1) = ',elec_num_tab(2) * (elec_num_tab(2) - 1) print*,'' print*,'ab_norm(istate) = ',ab_norm(istate) - print*,'N_alpha * N_beta = ',elec_num_tab(1) * elec_num_tab(2) + print*,'N_alpha * N_beta *2 = ',elec_num_tab(1) * elec_num_tab(2) * 2 print*,'' print*,'tot_norm(istate) = ',tot_norm(istate) - print*,'N(N-1)/2 = ',elec_num*(elec_num - 1)/2 + print*,'N(N-1)/2 = ',elec_num*(elec_num - 1) enddo + return wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 wee_ab_st_av = 0.d0 @@ -341,10 +342,10 @@ subroutine routine_full_mos rdm_ab_st_av = state_av_full_occ_2_rdm_ab_mo(l,k,j,i) rdm_tot_st_av = state_av_full_occ_2_rdm_spin_trace_mo(l,k,j,i) - wee_aa_st_av += vijkl * rdm_aa_st_av - wee_bb_st_av += vijkl * rdm_bb_st_av - wee_ab_st_av += vijkl * rdm_ab_st_av - wee_tot_st_av += vijkl * rdm_tot_st_av + wee_aa_st_av += 0.5d0 * vijkl * rdm_aa_st_av + wee_bb_st_av += 0.5d0 * vijkl * rdm_bb_st_av + wee_ab_st_av += 0.5d0 * vijkl * rdm_ab_st_av + wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av enddo enddo enddo diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index db793047..9e3d1df0 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -34,6 +34,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_ab_mo',wall_2 - wall_1 +! state_av_act_2_rdm_ab_mo *= 2.d0 END_PROVIDER @@ -48,7 +49,7 @@ ! ! 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_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC @@ -63,6 +64,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_aa_mo',wall_2 - wall_1 +! state_av_act_2_rdm_aa_mo *= 2.d0 END_PROVIDER @@ -76,7 +78,7 @@ ! ! 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_{\beta}^{act} * (N_{\beta}^{act} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC @@ -91,6 +93,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_bb_mo',wall_2 - wall_1 +! state_av_act_2_rdm_bb_mo *= 2.d0 END_PROVIDER @@ -104,7 +107,7 @@ ! ! 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} * (N_{elec} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index a4dea15f..7e35fc7b 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -29,7 +29,8 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)] enddo enddo enddo - two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) * 2.d0 + two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) +! * 2.d0 END_PROVIDER From fd63ab1355d481e90a3cef7401e934da63f04a19 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 17:33:43 +0100 Subject: [PATCH 12/16] state_av_full_occ_2_rdm_aa_mo work --- src/two_body_rdm/example.irp.f | 7 +- src/two_body_rdm/state_av_act_2rdm.irp.f | 6 +- .../state_av_full_orb_2_rdm.irp.f | 132 +++++++++--------- 3 files changed, 70 insertions(+), 75 deletions(-) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 985f6a5d..01e971ba 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -322,7 +322,7 @@ subroutine routine_full_mos print*,'N(N-1)/2 = ',elec_num*(elec_num - 1) enddo - return +! return wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 wee_ab_st_av = 0.d0 @@ -355,17 +355,12 @@ subroutine routine_full_mos print*,'' print*,'STATE AVERAGE ENERGY ' print*,'wee_aa_st_av = ',wee_aa_st_av - print*,'wee_aa_st_av_2 = ',wee_aa_st_av_2 print*,'wee_bb_st_av = ',wee_bb_st_av - print*,'wee_bb_st_av_2 = ',wee_bb_st_av_2 print*,'wee_ab_st_av = ',wee_ab_st_av - print*,'wee_ab_st_av_2 = ',wee_ab_st_av_2 print*,'Sum of components = ',wee_aa_st_av + wee_bb_st_av + wee_ab_st_av - print*,'Sum of components_2 = ',wee_aa_st_av_2 + wee_bb_st_av_2 + wee_ab_st_av_2 print*,'' print*,'Full energy ' print*,'wee_tot_st_av = ',wee_tot_st_av - print*,'wee_tot_st_av_2 = ',wee_tot_st_av_2 print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3 end diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index 9e3d1df0..a97ec3f9 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -34,7 +34,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_ab_mo',wall_2 - wall_1 -! state_av_act_2_rdm_ab_mo *= 2.d0 + state_av_act_2_rdm_ab_mo *= 2.d0 END_PROVIDER @@ -64,7 +64,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_aa_mo',wall_2 - wall_1 -! state_av_act_2_rdm_aa_mo *= 2.d0 + state_av_act_2_rdm_aa_mo *= 2.d0 END_PROVIDER @@ -93,7 +93,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_bb_mo',wall_2 - wall_1 -! state_av_act_2_rdm_bb_mo *= 2.d0 + state_av_act_2_rdm_bb_mo *= 2.d0 END_PROVIDER diff --git a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f index b3a5fe65..251b4d96 100644 --- a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f @@ -47,7 +47,7 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb) enddo enddo enddo @@ -61,7 +61,7 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -73,7 +73,7 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0 + state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0 enddo enddo @@ -90,7 +90,7 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb) enddo enddo enddo @@ -104,7 +104,7 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -116,7 +116,7 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0 + state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0 enddo enddo endif @@ -167,11 +167,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -181,8 +181,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo @@ -198,11 +198,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -212,8 +212,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo endif @@ -263,11 +263,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) enddo enddo enddo @@ -277,8 +277,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo @@ -294,11 +294,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) enddo enddo enddo @@ -308,8 +308,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo endif @@ -364,11 +364,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) enddo enddo enddo @@ -377,8 +377,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo if (.not.no_core_density)then @@ -390,11 +390,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) enddo enddo enddo @@ -403,8 +403,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo endif @@ -420,11 +420,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -433,8 +433,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo if (.not.no_core_density)then @@ -446,11 +446,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) ! 1 2 1 2 : EXCHANGE TERM - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -459,8 +459,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0 enddo enddo endif @@ -476,14 +476,14 @@ korb = list_inact(k) ! ALPHA INACTIVE - BETA ACTIVE ! alph beta alph beta - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) ! beta alph beta alph - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb) ! BETA INACTIVE - ALPHA ACTIVE ! beta alph beta alpha - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) ! alph beta alph beta - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -493,8 +493,8 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0 - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0 + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0d0 enddo enddo @@ -510,14 +510,14 @@ korb = list_core(k) !! BETA ACTIVE - ALPHA CORE ! alph beta alph beta - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb) ! beta alph beta alph - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb) !! ALPHA ACTIVE - BETA CORE ! alph beta alph beta - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb) ! beta alph beta alph - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb) enddo enddo enddo @@ -527,8 +527,8 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0 - state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0 + state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0D0 + state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0D0 enddo enddo From c95f1ee0ac22e79aa0f7b785e7503b1853079c80 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 17:42:44 +0100 Subject: [PATCH 13/16] changed all two-rdm with the normalization convtion to N(N-1) and not N(N-1)/2 --- src/two_body_rdm/example.irp.f | 4 +- src/two_body_rdm/full_orb_2_rdm.irp.f | 132 +++++++++++++------------- 2 files changed, 68 insertions(+), 68 deletions(-) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 01e971ba..30e2685a 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -150,7 +150,7 @@ subroutine routine_active_only print*,'accu_ab = ',accu_ab print*,'N_a N_b = ', elec_beta_num*elec_alpha_num print*,'accu_tot = ',accu_tot - print*,'Ne(Ne-1)/2 = ',(elec_num-1)*elec_num + print*,'Ne(Ne-1) = ',(elec_num-1)*elec_num enddo wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 @@ -319,7 +319,7 @@ subroutine routine_full_mos print*,'N_alpha * N_beta *2 = ',elec_num_tab(1) * elec_num_tab(2) * 2 print*,'' print*,'tot_norm(istate) = ',tot_norm(istate) - print*,'N(N-1)/2 = ',elec_num*(elec_num - 1) + print*,'N(N-1) = ',elec_num*(elec_num - 1) enddo ! return diff --git a/src/two_body_rdm/full_orb_2_rdm.irp.f b/src/two_body_rdm/full_orb_2_rdm.irp.f index fba88172..a3de8ea9 100644 --- a/src/two_body_rdm/full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/full_orb_2_rdm.irp.f @@ -50,7 +50,7 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -64,7 +64,7 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -76,7 +76,7 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0 + full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0 enddo enddo @@ -93,7 +93,7 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -107,7 +107,7 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -119,7 +119,7 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0 + full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0 enddo enddo endif @@ -172,11 +172,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -186,8 +186,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo @@ -203,11 +203,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -217,8 +217,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo endif @@ -270,11 +270,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -284,8 +284,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo @@ -301,11 +301,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -315,8 +315,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo endif @@ -377,11 +377,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -390,8 +390,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo if (.not.no_core_density)then @@ -403,11 +403,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -416,8 +416,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo endif @@ -433,11 +433,11 @@ do k = 1, n_inact_orb korb = list_inact(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -446,8 +446,8 @@ jorb = list_inact(j) do k = 1, n_inact_orb korb = list_inact(k) - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo if (.not.no_core_density)then @@ -459,11 +459,11 @@ do k = 1, n_core_orb korb = list_core(k) ! 1 2 1 2 : DIRECT TERM - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) ! 1 2 1 2 : EXCHANGE TERM - full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -472,8 +472,8 @@ jorb = list_core(j) do k = 1, n_core_orb korb = list_core(k) - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0 enddo enddo endif @@ -489,14 +489,14 @@ korb = list_inact(k) ! ALPHA INACTIVE - BETA ACTIVE ! alph beta alph beta - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) ! beta alph beta alph - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate) ! BETA INACTIVE - ALPHA ACTIVE ! beta alph beta alpha - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) ! alph beta alph beta - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -506,8 +506,8 @@ do k = 1, n_inact_orb korb = list_inact(k) ! alph beta alph beta - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0 - full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0 + full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0 enddo enddo @@ -523,14 +523,14 @@ korb = list_core(k) !! BETA ACTIVE - ALPHA CORE ! alph beta alph beta - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate) ! beta alph beta alph - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate) !! ALPHA ACTIVE - BETA CORE ! alph beta alph beta - full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate) ! beta alph beta alph - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -540,8 +540,8 @@ do k = 1, n_core_orb korb = list_core(k) ! alph beta alph beta - full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0 - full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0 + full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0 + full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0 enddo enddo From 8515fcf93fae51d4658f2a9911f39d6bf1a2fced Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 17:48:42 +0100 Subject: [PATCH 14/16] updated DOC for correct normalization of two-rdms --- src/two_body_rdm/act_2_rdm.irp.f | 19 +++++++--------- src/two_body_rdm/full_orb_2_rdm.irp.f | 22 ++++++++----------- src/two_body_rdm/state_av_act_2rdm.irp.f | 12 +++++----- .../state_av_full_orb_2_rdm.irp.f | 16 +++++--------- 4 files changed, 27 insertions(+), 42 deletions(-) diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 61ae4e47..c550e991 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -4,21 +4,18 @@ BEGIN_DOC ! 12 12 ! 1 2 1 2 == -! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons +! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons ! -! +! +! +! + ! ! 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_{\alpha}^{act} * N_{\beta}^{act} +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2 ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! -! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta -! -! act_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta -! -! Therefore you don't necessary have symmetry between electron 1 and 2 END_DOC integer :: ispin double precision :: wall_1, wall_2 @@ -57,7 +54,7 @@ ! ! 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_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC @@ -98,7 +95,7 @@ ! ! 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_{\beta}^{act} * (N_{\beta}^{act} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC @@ -138,7 +135,7 @@ ! ! 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)/2 +! 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 diff --git a/src/two_body_rdm/full_orb_2_rdm.irp.f b/src/two_body_rdm/full_orb_2_rdm.irp.f index a3de8ea9..78a28acf 100644 --- a/src/two_body_rdm/full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/full_orb_2_rdm.irp.f @@ -4,22 +4,18 @@ full_occ_2_rdm_ab_mo = 0.d0 integer :: i,j,k,l,iorb,jorb,korb,lorb,istate BEGIN_DOC -! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons +! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta + beta/alpha electrons ! -! +! +! +! + ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2 ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! -! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA -! -! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta -! -! Therefore you don't necessary have symmetry between electron 1 and 2 -! ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO ARE SET TO ZERO END_DOC full_occ_2_rdm_ab_mo = 0.d0 @@ -139,7 +135,7 @@ ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! @@ -237,7 +233,7 @@ ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! @@ -335,14 +331,14 @@ ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero ! The two-electron energy of each state can be computed as: ! -! \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll > +! \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * 1/2 * < ii jj | kk ll > ! ! with ii = list_core_inact_act(i), jj = list_core_inact_act(j), kk = list_core_inact_act(k), ll = list_core_inact_act(l) END_DOC diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index a97ec3f9..cd417a9d 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -2,21 +2,16 @@ implicit none double precision, allocatable :: state_weights(:) BEGIN_DOC -! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs +! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha/beta+beta/alpha electron pairs ! ! = \sum_{istate} w(istate) * ! ! 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_{\alpha}^{act} * N_{\beta}^{act} +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2 ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! -! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta -! -! state_av_act_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta -! -! Therefore you don't necessary have symmetry between electron 1 and 2 END_DOC allocate(state_weights(N_states)) state_weights = state_average_weight @@ -34,6 +29,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_ab_mo',wall_2 - wall_1 + ! factor 2 to have the correct normalization factor state_av_act_2_rdm_ab_mo *= 2.d0 END_PROVIDER @@ -64,6 +60,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_aa_mo',wall_2 - wall_1 + ! factor 2 to have the correct normalization factor state_av_act_2_rdm_aa_mo *= 2.d0 END_PROVIDER @@ -93,6 +90,7 @@ call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Wall time to provide state_av_act_2_rdm_bb_mo',wall_2 - wall_1 + ! factor 2 to have the correct normalization factor state_av_act_2_rdm_bb_mo *= 2.d0 END_PROVIDER diff --git a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f index 251b4d96..2e44665d 100644 --- a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f @@ -4,22 +4,16 @@ state_av_full_occ_2_rdm_ab_mo = 0.d0 integer :: i,j,k,l,iorb,jorb,korb,lorb BEGIN_DOC -! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta electrons +! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta + beta/alpha electrons ! ! = \sum_{istate} w(istate) * ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2 ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! -! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA -! -! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta -! -! Therefore you don't necessary have symmetry between electron 1 and 2 -! ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero END_DOC state_av_full_occ_2_rdm_ab_mo = 0.d0 @@ -135,7 +129,7 @@ ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! @@ -231,7 +225,7 @@ ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! @@ -328,7 +322,7 @@ ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active ! -! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2 +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1) ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! From b51d72d5b905ef5ac198d42cc4b53ca90192d90d Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 18:35:51 +0100 Subject: [PATCH 15/16] changed the factor 2 in basis_correction and mu_of_r in order to adapt to new normalization factor --- src/basis_correction/pbe_on_top.irp.f | 24 +++++++++++++++--------- src/mu_of_r/f_psi_i_a_v_utils.irp.f | 2 +- src/mu_of_r/f_psi_utils.irp.f | 2 +- src/mu_of_r/mu_of_r_conditions.irp.f | 2 +- 4 files changed, 18 insertions(+), 12 deletions(-) diff --git a/src/basis_correction/pbe_on_top.irp.f b/src/basis_correction/pbe_on_top.irp.f index cc9cdec9..a25fd61b 100644 --- a/src/basis_correction/pbe_on_top.irp.f +++ b/src/basis_correction/pbe_on_top.irp.f @@ -41,13 +41,15 @@ if(mu_of_r_potential == "cas_ful")then ! You take the on-top of the CAS wave function which is computed with mu(r) - on_top = on_top_cas_mu_r(ipoint,istate) + ! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r + on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate) else ! You take the on-top of the CAS wave function computed separately + ! No factor 2 because convention N(N-1) in provider total_cas_on_top_density on_top = total_cas_on_top_density(ipoint,istate) endif -! We take the extrapolated on-top pair density * 2 because of normalization - on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top) +! We take the extrapolated on-top pair density + on_top_extrap = mu_correction_of_on_top(mu,on_top) call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE) @@ -103,13 +105,15 @@ if(mu_of_r_potential == "cas_ful")then ! You take the on-top of the CAS wave function which is computed with mu(r) - on_top = on_top_cas_mu_r(ipoint,istate) + ! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r + on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate) else ! You take the on-top of the CAS wave function computed separately + ! No factor 2 because convention N(N-1) in provider total_cas_on_top_density on_top = total_cas_on_top_density(ipoint,istate) endif -! We take the extrapolated on-top pair density * 2 because of normalization - on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top) +! We take the extrapolated on-top pair density + on_top_extrap = mu_correction_of_on_top(mu,on_top) call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE) @@ -165,13 +169,15 @@ if(mu_of_r_potential == "cas_ful")then ! You take the on-top of the CAS wave function which is computed with mu(r) - on_top = on_top_cas_mu_r(ipoint,istate) + ! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r + on_top = 1.d0 * on_top_cas_mu_r(ipoint,istate) else ! You take the on-top of the CAS wave function computed separately + ! No factor 2 because convention N(N-1) in provider total_cas_on_top_density on_top = total_cas_on_top_density(ipoint,istate) endif -! We DO NOT take the extrapolated on-top pair density, but there is * 2 because of normalization - on_top_extrap = 2.d0 * on_top +! We DO NOT take the extrapolated on-top pair density + on_top_extrap = on_top call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE) diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f index 427da199..0d08e193 100644 --- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f +++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f @@ -194,7 +194,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) do b = 1, n_act_orb ! 2 do c = 1, n_act_orb ! 1 do d = 1, n_act_orb ! 2 - rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) + rho = mos_array_act_r1(c) * mos_array_act_r2(d) * 0.5d0 * act_2_rdm_ab_mo(d,c,b,a,istate) rho_tilde(b,a) += rho two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) enddo diff --git a/src/mu_of_r/f_psi_utils.irp.f b/src/mu_of_r/f_psi_utils.irp.f index bdd76f18..95f10f36 100644 --- a/src/mu_of_r/f_psi_utils.irp.f +++ b/src/mu_of_r/f_psi_utils.irp.f @@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd_trans, (n_points_final_ do j = 1, n_core_inact_act_orb do i = 1, n_core_inact_act_orb ! 1 2 1 2 - full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint) + full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += 0.5d0 * full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint) enddo enddo enddo diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index 5c41acdc..f9c3b3b3 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -110,7 +110,7 @@ do istate = 1, N_states do ipoint = 1, n_points_final_grid f_psi = f_psi_cas_ab(ipoint,istate) - on_top = on_top_cas_mu_r(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then w_psi = 1.d+10 else From 46dfb551a8562feb693c2c1d6190ba63be779f23 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 4 Mar 2023 13:07:12 +0100 Subject: [PATCH 16/16] modified the factor two in rdm --- src/basis_correction/pbe_on_top.irp.f | 12 +++--------- src/mu_of_r/f_hf_utils.irp.f | 5 +++++ src/mu_of_r/f_psi_i_a_v_utils.irp.f | 10 +++++++++- 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/basis_correction/pbe_on_top.irp.f b/src/basis_correction/pbe_on_top.irp.f index a25fd61b..9167f459 100644 --- a/src/basis_correction/pbe_on_top.irp.f +++ b/src/basis_correction/pbe_on_top.irp.f @@ -41,11 +41,9 @@ if(mu_of_r_potential == "cas_ful")then ! You take the on-top of the CAS wave function which is computed with mu(r) - ! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r - on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) else ! You take the on-top of the CAS wave function computed separately - ! No factor 2 because convention N(N-1) in provider total_cas_on_top_density on_top = total_cas_on_top_density(ipoint,istate) endif ! We take the extrapolated on-top pair density @@ -105,11 +103,9 @@ if(mu_of_r_potential == "cas_ful")then ! You take the on-top of the CAS wave function which is computed with mu(r) - ! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r - on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) else ! You take the on-top of the CAS wave function computed separately - ! No factor 2 because convention N(N-1) in provider total_cas_on_top_density on_top = total_cas_on_top_density(ipoint,istate) endif ! We take the extrapolated on-top pair density @@ -169,11 +165,9 @@ if(mu_of_r_potential == "cas_ful")then ! You take the on-top of the CAS wave function which is computed with mu(r) - ! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r - on_top = 1.d0 * on_top_cas_mu_r(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) else ! You take the on-top of the CAS wave function computed separately - ! No factor 2 because convention N(N-1) in provider total_cas_on_top_density on_top = total_cas_on_top_density(ipoint,istate) endif ! We DO NOT take the extrapolated on-top pair density diff --git a/src/mu_of_r/f_hf_utils.irp.f b/src/mu_of_r/f_hf_utils.irp.f index 8480a288..102b40a0 100644 --- a/src/mu_of_r/f_hf_utils.irp.f +++ b/src/mu_of_r/f_hf_utils.irp.f @@ -86,6 +86,9 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens) enddo enddo enddo + ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm + f_HF_val_ab *= 2.d0 + two_bod_dens *= 2.d0 end @@ -136,4 +139,6 @@ subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab) enddo enddo enddo + ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm + int_f_HF_val_ab *= 2.d0 end diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f index 0d08e193..69fa16ff 100644 --- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f +++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f @@ -49,6 +49,9 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens) enddo enddo enddo + ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm + f_ii_val_ab *= 2.d0 + two_bod_dens *= 2.d0 end @@ -142,6 +145,9 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) f_ia_val_ab += v_tilde(i,a) * rho_tilde(i,a) enddo enddo + ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm + f_ia_val_ab *= 2.d0 + two_bod_dens *= 2.d0 end @@ -194,7 +200,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) do b = 1, n_act_orb ! 2 do c = 1, n_act_orb ! 1 do d = 1, n_act_orb ! 2 - rho = mos_array_act_r1(c) * mos_array_act_r2(d) * 0.5d0 * act_2_rdm_ab_mo(d,c,b,a,istate) + rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) rho_tilde(b,a) += rho two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) enddo @@ -222,6 +228,8 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) f_aa_val_ab += v_tilde(b,a) * rho_tilde(b,a) enddo enddo + + ! DO NOT multiply by two as in give_f_ii_val_ab and give_f_ia_val_ab because the N(N-1) normalization condition of the active two-rdm end BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)]