LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

Merge pull request #11 from pfloos/work_env

Browse Source 
Work env
This commit is contained in:
AbdAmmar 2024-11-15 11:38:01 +01:00 committed by GitHub
commit 1cccc81c0d
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
10 changed files with 685 additions and 476 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
PyDuck.py
View File

@ -1,6 +1,6 @@
#!/usr/bin/env python #!/usr/bin/env python
import os import os, sys
import argparse import argparse
import pyscf import pyscf
from pyscf import gto from pyscf import gto
@ -8,6 +8,10 @@ import numpy as np
import subprocess import subprocess
#Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository #Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository
if "QUACK_ROOT" not in os.environ:
print("Please set the QUACK_ROOT environment variable, for example:\n")
print("$ export QUACK_ROOT={0}".format(os.getcwd()))
sys.exit(1)
QuAcK_dir=os.environ.get('QUACK_ROOT','./') QuAcK_dir=os.environ.get('QUACK_ROOT','./')
#Create the argument parser object and gives a description of the script #Create the argument parser object and gives a description of the script
@ -37,7 +41,7 @@ print_2e=args.print_2e
working_dir=args.working_dir working_dir=args.working_dir
#Read molecule #Read molecule
f = open(QuAcK_dir + '/mol/' + xyz,'r') f = open(working_dir+'/mol/'+xyz,'r')
lines = f.read().splitlines() lines = f.read().splitlines()
nbAt = int(lines.pop(0)) nbAt = int(lines.pop(0))
lines.pop(0) lines.pop(0)
@ -70,6 +74,7 @@ nelec=mol.nelec #Access the number of electrons
nalpha=nelec[0] nalpha=nelec[0]
nbeta=nelec[1] nbeta=nelec[1]
subprocess.call(['mkdir', '-p', working_dir+'/input'])
f = open(working_dir+'/input/molecule','w') f = open(working_dir+'/input/molecule','w')
f.write('# nAt nEla nElb nCore nRyd\n') f.write('# nAt nEla nElb nCore nRyd\n')
f.write(str(mol.natm)+' '+str(nalpha)+' '+str(nbeta)+' '+str(0)+' '+str(0)+'\n') f.write(str(mol.natm)+' '+str(nalpha)+' '+str(nbeta)+' '+str(0)+' '+str(0)+'\n')
@ -79,7 +84,8 @@ for i in range(len(list_pos_atom)):
f.close() f.close()
#Compute nuclear energy and put it in a file #Compute nuclear energy and put it in a file
subprocess.call(['rm', working_dir + '/int/ENuc.dat']) subprocess.call(['mkdir', '-p', working_dir+'/int'])
subprocess.call(['rm', '-f', working_dir + '/int/ENuc.dat'])
f = open(working_dir+'/int/ENuc.dat','w') f = open(working_dir+'/int/ENuc.dat','w')
f.write(str(mol.energy_nuc())) f.write(str(mol.energy_nuc()))
f.write(' ') f.write(' ')
@ -93,14 +99,14 @@ dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators
x,y,z = dipole[0],dipole[1],dipole[2] x,y,z = dipole[0],dipole[1],dipole[2]
norb = len(ovlp) # nBAS_AOs norb = len(ovlp) # nBAS_AOs
subprocess.call(['rm', working_dir + '/int/nBas.dat']) subprocess.call(['rm', '-f', working_dir + '/int/nBas.dat'])
f = open(working_dir+'/int/nBas.dat','w') f = open(working_dir+'/int/nBas.dat','w')
f.write(" {} ".format(str(norb))) f.write(" {} ".format(str(norb)))
f.close() f.close()
def write_matrix_to_file(matrix,size,file,cutoff=1e-15): def write_matrix_to_file(matrix,size,file,cutoff=1e-15):
f = open(file, 'a') f = open(file, 'w')
for i in range(size): for i in range(size):
for j in range(i,size): for j in range(i,size):
if abs(matrix[i][j]) > cutoff: if abs(matrix[i][j]) > cutoff:
@ -110,23 +116,23 @@ def write_matrix_to_file(matrix,size,file,cutoff=1e-15):
#Write all 1 electron quantities in files #Write all 1 electron quantities in files
#Ov,Nuc,Kin,x,y,z #Ov,Nuc,Kin,x,y,z
subprocess.call(['rm', working_dir + '/int/Ov.dat']) subprocess.call(['rm', '-f', working_dir + '/int/Ov.dat'])
write_matrix_to_file(ovlp,norb,working_dir+'/int/Ov.dat') write_matrix_to_file(ovlp,norb,working_dir+'/int/Ov.dat')
subprocess.call(['rm', working_dir + '/int/Nuc.dat']) subprocess.call(['rm', '-f', working_dir + '/int/Nuc.dat'])
write_matrix_to_file(v1e,norb,working_dir+'/int/Nuc.dat') write_matrix_to_file(v1e,norb,working_dir+'/int/Nuc.dat')
subprocess.call(['rm', working_dir + '/int/Kin.dat']) subprocess.call(['rm', '-f', working_dir + '/int/Kin.dat'])
write_matrix_to_file(t1e,norb,working_dir+'/int/Kin.dat') write_matrix_to_file(t1e,norb,working_dir+'/int/Kin.dat')
subprocess.call(['rm', working_dir + '/int/x.dat']) subprocess.call(['rm', '-f', working_dir + '/int/x.dat'])
write_matrix_to_file(x,norb,working_dir+'/int/x.dat') write_matrix_to_file(x,norb,working_dir+'/int/x.dat')
subprocess.call(['rm', working_dir + '/int/y.dat']) subprocess.call(['rm', '-f', working_dir + '/int/y.dat'])
write_matrix_to_file(y,norb,working_dir+'/int/y.dat') write_matrix_to_file(y,norb,working_dir+'/int/y.dat')
subprocess.call(['rm', working_dir + '/int/z.dat']) subprocess.call(['rm', '-f', working_dir + '/int/z.dat'])
write_matrix_to_file(z,norb,working_dir+'/int/z.dat') write_matrix_to_file(z,norb,working_dir+'/int/z.dat')
eri_ao = mol.intor('int2e') eri_ao = mol.intor('int2e')
def write_tensor_to_file(tensor,size,file,cutoff=1e-15): def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
f = open(file, 'a') f = open(file, 'w')
for i in range(size): for i in range(size):
for j in range(i,size): for j in range(i,size):
for k in range(i,size): for k in range(i,size):
@ -140,15 +146,18 @@ def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
# Write two-electron integrals # Write two-electron integrals
if print_2e: if print_2e:
# (formatted) # (formatted)
subprocess.call(['rm', working_dir + '/int/ERI.dat']) subprocess.call(['rm', '-f', working_dir + '/int/ERI.dat'])
write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat') write_tensor_to_file(eri_ao, norb, working_dir + '/int/ERI.dat')
else: else:
# (binary) # (binary)
subprocess.call(['rm', working_dir + '/int/ERI.bin']) subprocess.call(['rm', '-f', working_dir + '/int/ERI.bin'])
# chem -> phys notation # chem -> phys notation
eri_ao = eri_ao.transpose(0, 2, 1, 3) eri_ao = eri_ao.transpose(0, 2, 1, 3)
eri_ao.tofile('int/ERI.bin') f = open(working_dir + '/int/ERI.bin', 'w')
eri_ao.tofile(working_dir + '/int/ERI.bin')
f.close()
#Execute the QuAcK fortran program #Execute the QuAcK fortran program
subprocess.call(QuAcK_dir+'/bin/QuAcK') subprocess.call([QuAcK_dir + '/bin/QuAcK', working_dir])

26
src/QuAcK/QuAcK.f90
View File

@ -71,6 +71,16 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest logical :: dotest,doRtest,doUtest,doGtest
character(len=256) :: working_dir
! Check if the right number of arguments is provided
if(command_argument_count() < 1) then
print *, "No working directory provided."
stop
else
call get_command_argument(1, working_dir)
endif
!-------------! !-------------!
! Hello World ! ! Hello World !
!-------------! !-------------!
@ -95,7 +105,8 @@ program QuAcK
! Method selection ! ! Method selection !
!------------------! !------------------!
call read_methods(doRHF,doUHF,doGHF,doROHF, & call read_methods(working_dir, &
doRHF,doUHF,doGHF,doROHF, &
doMP2,doMP3, & doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD, & dodrCCD,dorCCD,docrCCD,dolCCD, &
@ -112,7 +123,8 @@ program QuAcK
! Read options for methods ! ! Read options for methods !
!--------------------------! !--------------------------!
call read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, & call read_options(working_dir, &
maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
reg_MP, & reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, & maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, & TDA,spin_conserved,spin_flip, &
@ -133,18 +145,18 @@ program QuAcK
! nOrb = number of orbitals ! ! nOrb = number of orbitals !
!------------------------------------! !------------------------------------!
call read_molecule(nNuc,nO,nC,nR) call read_molecule(working_dir,nNuc,nO,nC,nR)
allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
! Read geometry ! Read geometry
call read_geometry(nNuc,ZNuc,rNuc,ENuc) call read_geometry(working_dir,nNuc,ZNuc,rNuc,ENuc)
!---------------------------------------! !---------------------------------------!
! Read basis set information from PySCF ! ! Read basis set information from PySCF !
!---------------------------------------! !---------------------------------------!
call read_basis_pyscf(nBas,nO,nV) call read_basis_pyscf(working_dir,nBas,nO,nV)
!--------------------------------------! !--------------------------------------!
! Read one- and two-electron integrals ! ! Read one- and two-electron integrals !
@ -163,8 +175,8 @@ program QuAcK
call wall_time(start_int) call wall_time(start_int)
call read_integrals(nBas,S,T,V,Hc,ERI_AO) call read_integrals(working_dir,nBas,S,T,V,Hc,ERI_AO)
call read_dipole_integrals(nBas,dipole_int_AO) call read_dipole_integrals(working_dir,nBas,dipole_int_AO)
call wall_time(end_int) call wall_time(end_int)

271
src/QuAcK/read_methods.f90
View File

@ -1,4 +1,5 @@
subroutine read_methods(doRHF,doUHF,doGHF,doROHF, & subroutine read_methods(working_dir, &
doRHF,doUHF,doGHF,doROHF, &
doMP2,doMP3, & doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
@ -17,6 +18,10 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF, &
! Input variables ! Input variables
character(len=256),intent(in) :: working_dir
! Output variables
logical,intent(out) :: doRHF,doUHF,doGHF,doROHF logical,intent(out) :: doRHF,doUHF,doGHF,doROHF
logical,intent(out) :: doMP2,doMP3 logical,intent(out) :: doMP2,doMP3
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
@ -33,170 +38,180 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF, &
! Local variables ! Local variables
character(len=1) :: ans1,ans2,ans3,ans4,ans5,ans6 character(len=1) :: ans1,ans2,ans3,ans4,ans5,ans6
integer :: status
character(len=256) :: file_path
! Open file with method specification
open(unit=1,file='input/methods') file_path = trim(working_dir) // '/input/methods'
open(unit=1, file=file_path, status='old', action='read', iostat=status)
! Read mean-field methods if(status /= 0) then
doRHF = .false. print *, "Error opening file: ", file_path
doUHF = .false. stop
doGHF = .false.
doROHF = .false.
read(1,*) else
read(1,*) ans1,ans2,ans3,ans4
if(ans1 == 'T') doRHF = .true.
if(ans2 == 'T') doUHF = .true.
if(ans3 == 'T') doGHF = .true.
if(ans4 == 'T') doROHF = .true.
! Read MPn methods ! Read mean-field methods
doMP2 = .false. doRHF = .false.
doMP3 = .false. doUHF = .false.
doGHF = .false.
doROHF = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2 read(1,*) ans1,ans2,ans3,ans4
if(ans1 == 'T') doMP2 = .true. if(ans1 == 'T') doRHF = .true.
if(ans2 == 'T') doMP3 = .true. if(ans2 == 'T') doUHF = .true.
if(ans3 == 'T') doGHF = .true.
if(ans4 == 'T') doROHF = .true.
! Read CC methods ! Read MPn methods
doCCD = .false. doMP2 = .false.
dopCCD = .false. doMP3 = .false.
doDCD = .false.
doCCSD = .false.
doCCSDT = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4,ans5 read(1,*) ans1,ans2
if(ans1 == 'T') doCCD = .true. if(ans1 == 'T') doMP2 = .true.
if(ans2 == 'T') dopCCD = .true. if(ans2 == 'T') doMP3 = .true.
if(ans3 == 'T') doDCD = .true.
if(ans4 == 'T') doCCSD = .true.
if(ans5 == 'T') doCCSDT = .true.
! Read weird CC methods ! Read CC methods
do_drCCD = .false. doCCD = .false.
do_rCCD = .false. dopCCD = .false.
do_crCCD = .false. doDCD = .false.
do_lCCD = .false. doCCSD = .false.
doCCSDT = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4 read(1,*) ans1,ans2,ans3,ans4,ans5
if(ans1 == 'T') do_drCCD = .true. if(ans1 == 'T') doCCD = .true.
if(ans2 == 'T') do_rCCD = .true. if(ans2 == 'T') dopCCD = .true.
if(ans3 == 'T') do_crCCD = .true. if(ans3 == 'T') doDCD = .true.
if(ans4 == 'T') do_lCCD = .true. if(ans4 == 'T') doCCSD = .true.
if(ans5 == 'T') doCCSDT = .true.
! Read CI methods ! Read weird CC methods
doCIS = .false. do_drCCD = .false.
doCIS_D = .false. do_rCCD = .false.
doCID = .false. do_crCCD = .false.
doCISD = .false. do_lCCD = .false.
doFCI = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4,ans5 read(1,*) ans1,ans2,ans3,ans4
if(ans1 == 'T') doCIS = .true. if(ans1 == 'T') do_drCCD = .true.
if(ans2 == 'T') doCIS_D = .true. if(ans2 == 'T') do_rCCD = .true.
if(ans3 == 'T') doCID = .true. if(ans3 == 'T') do_crCCD = .true.
if(ans4 == 'T') doCISD = .true. if(ans4 == 'T') do_lCCD = .true.
if(ans5 == 'T') doFCI = .true.
if(doCIS_D) doCIS = .true.
! Read RPA methods ! Read CI methods
dophRPA = .false. doCIS = .false.
dophRPAx = .false. doCIS_D = .false.
docrRPA = .false. doCID = .false.
doppRPA = .false. doCISD = .false.
doFCI = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4 read(1,*) ans1,ans2,ans3,ans4,ans5
if(ans1 == 'T') dophRPA = .true. if(ans1 == 'T') doCIS = .true.
if(ans2 == 'T') dophRPAx = .true. if(ans2 == 'T') doCIS_D = .true.
if(ans3 == 'T') docrRPA = .true. if(ans3 == 'T') doCID = .true.
if(ans4 == 'T') doppRPA = .true. if(ans4 == 'T') doCISD = .true.
if(ans5 == 'T') doFCI = .true.
if(doCIS_D) doCIS = .true.
! Read Green's function methods ! Read RPA methods
doG0F2 = .false. dophRPA = .false.
doevGF2 = .false. dophRPAx = .false.
doqsGF2 = .false. docrRPA = .false.
doufG0F02 = .false. doppRPA = .false.
doG0F3 = .false.
doevGF3 = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4,ans5,ans6 read(1,*) ans1,ans2,ans3,ans4
if(ans1 == 'T') doG0F2 = .true. if(ans1 == 'T') dophRPA = .true.
if(ans2 == 'T') doevGF2 = .true. if(ans2 == 'T') dophRPAx = .true.
if(ans3 == 'T') doqsGF2 = .true. if(ans3 == 'T') docrRPA = .true.
if(ans4 == 'T') doufG0F02 = .true. if(ans4 == 'T') doppRPA = .true.
if(ans5 == 'T') doG0F3 = .true.
if(ans6 == 'T') doevGF3 = .true.
! Read GW methods ! Read Green's function methods
doG0W0 = .false. doG0F2 = .false.
doevGW = .false. doevGF2 = .false.
doqsGW = .false. doqsGF2 = .false.
doufG0W0 = .false. doufG0F02 = .false.
doufGW = .false. doG0F3 = .false.
doevGF3 = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4,ans5 read(1,*) ans1,ans2,ans3,ans4,ans5,ans6
if(ans1 == 'T') doG0W0 = .true. if(ans1 == 'T') doG0F2 = .true.
if(ans2 == 'T') doevGW = .true. if(ans2 == 'T') doevGF2 = .true.
if(ans3 == 'T') doqsGW = .true. if(ans3 == 'T') doqsGF2 = .true.
if(ans4 == 'T') doufG0W0 = .true. if(ans4 == 'T') doufG0F02 = .true.
if(ans5 == 'T') doufGW = .true. if(ans5 == 'T') doG0F3 = .true.
if(ans6 == 'T') doevGF3 = .true.
! Read GTpp methods ! Read GW methods
doG0T0pp = .false. doG0W0 = .false.
doevGTpp = .false. doevGW = .false.
doqsGTpp = .false. doqsGW = .false.
doufG0T0pp = .false. doufG0W0 = .false.
doufGW = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4 read(1,*) ans1,ans2,ans3,ans4,ans5
if(ans1 == 'T') doG0T0pp = .true. if(ans1 == 'T') doG0W0 = .true.
if(ans2 == 'T') doevGTpp = .true. if(ans2 == 'T') doevGW = .true.
if(ans3 == 'T') doqsGTpp = .true. if(ans3 == 'T') doqsGW = .true.
if(ans4 == 'T') doufG0T0pp = .true. if(ans4 == 'T') doufG0W0 = .true.
if(ans5 == 'T') doufGW = .true.
! Read GTeh methods ! Read GTpp methods
doG0T0eh = .false. doG0T0pp = .false.
doevGTeh = .false. doevGTpp = .false.
doqsGTeh = .false. doqsGTpp = .false.
doufG0T0pp = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3 read(1,*) ans1,ans2,ans3,ans4
if(ans1 == 'T') doG0T0eh = .true. if(ans1 == 'T') doG0T0pp = .true.
if(ans2 == 'T') doevGTeh = .true. if(ans2 == 'T') doevGTpp = .true.
if(ans3 == 'T') doqsGTeh = .true. if(ans3 == 'T') doqsGTpp = .true.
if(ans4 == 'T') doufG0T0pp = .true.
! Read test ! Read GTeh methods
doRtest = .false. doG0T0eh = .false.
doUtest = .false. doevGTeh = .false.
doGtest = .false. doqsGTeh = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3 read(1,*) ans1,ans2,ans3
if(ans1 == 'T') doRtest = .true. if(ans1 == 'T') doG0T0eh = .true.
if(ans2 == 'T') doUtest = .true. if(ans2 == 'T') doevGTeh = .true.
if(ans3 == 'T') doGtest = .true. if(ans3 == 'T') doqsGTeh = .true.
! Close file ! Read test
doRtest = .false.
doUtest = .false.
doGtest = .false.
read(1,*)
read(1,*) ans1,ans2,ans3
if(ans1 == 'T') doRtest = .true.
if(ans2 == 'T') doUtest = .true.
if(ans3 == 'T') doGtest = .true.
endif
! Close file
close(unit=1) close(unit=1)
end subroutine end subroutine

208
src/QuAcK/read_options.f90
View File

@ -1,4 +1,5 @@
subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, & subroutine read_options(working_dir, &
maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
reg_MP, & reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, & maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, & TDA,spin_conserved,spin_flip, &
@ -14,6 +15,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
! Input variables ! Input variables
character(len=256),intent(in) :: working_dir
! Output variables
integer,intent(out) :: maxSCF_HF integer,intent(out) :: maxSCF_HF
double precision,intent(out) :: thresh_HF double precision,intent(out) :: thresh_HF
integer,intent(out) :: max_diis_HF integer,intent(out) :: max_diis_HF
@ -70,140 +75,151 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
! Local variables ! Local variables
character(len=1) :: ans1,ans2,ans3,ans4,ans5 character(len=1) :: ans1,ans2,ans3,ans4,ans5
integer :: status
character(len=256) :: file_path
! Open file with method specification ! Open file with method specification
open(unit=1,file='input/options') file_path = trim(working_dir) // '/input/options'
open(unit=1, file=file_path, status='old', action='read', iostat=status)
! Read HF options if(status /= 0) then
maxSCF_HF = 64 print *, "Error opening file: ", file_path
thresh_HF = 1d-6 stop
max_diis_HF = 1
guess_type = 1
mix = 0d0
level_shift = 0d0
dostab = .false.
dosearch = .false.
read(1,*) else
read(1,*) maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,ans1,ans2
if(ans1 == 'T') dostab = .true. ! Read HF options
if(ans2 == 'T') dosearch = .true.
! Read MPn options maxSCF_HF = 64
thresh_HF = 1d-6
max_diis_HF = 1
guess_type = 1
mix = 0d0
level_shift = 0d0
dostab = .false.
dosearch = .false.
reg_MP = .false. read(1,*)
read(1,*) read(1,*) maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,ans1,ans2
read(1,*) ans1
if(ans1 == 'T') reg_MP = .true. if(ans1 == 'T') dostab = .true.
if(ans2 == 'T') dosearch = .true.
! Read CC options ! Read MPn options
maxSCF_CC = 64 reg_MP = .false.
thresh_CC = 1d-5 read(1,*)
max_diis_CC = 1 read(1,*) ans1
read(1,*) if(ans1 == 'T') reg_MP = .true.
read(1,*) maxSCF_CC,thresh_CC,max_diis_CC
! Read excited state options ! Read CC options
TDA = .false. maxSCF_CC = 64
spin_conserved = .false. thresh_CC = 1d-5
spin_flip = .false. max_diis_CC = 1
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3 read(1,*) maxSCF_CC,thresh_CC,max_diis_CC
if(ans1 == 'T') TDA = .true. ! Read excited state options
if(ans2 == 'T') spin_conserved = .true.
if(ans3 == 'T') spin_flip = .true.
! Read GF options TDA = .false.
spin_conserved = .false.
spin_flip = .false.
maxSCF_GF = 64 read(1,*)
thresh_GF = 1d-5 read(1,*) ans1,ans2,ans3
max_diis_GF = 1
lin_GF = .false.
eta_GF = 0d0
renorm_GF = 0
reg_GF = .false.
read(1,*) if(ans1 == 'T') TDA = .true.
read(1,*) maxSCF_GF,thresh_GF,max_diis_GF,ans1,eta_GF,renorm_GF,ans2 if(ans2 == 'T') spin_conserved = .true.
if(ans3 == 'T') spin_flip = .true.
if(ans1 == 'T') lin_GF = .true. ! Read GF options
if(ans2 == 'T') reg_GF = .true.
! Read GW options maxSCF_GF = 64
thresh_GF = 1d-5
max_diis_GF = 1
lin_GF = .false.
eta_GF = 0d0
renorm_GF = 0
reg_GF = .false.
maxSCF_GW = 64 read(1,*)
thresh_GW = 1d-5 read(1,*) maxSCF_GF,thresh_GF,max_diis_GF,ans1,eta_GF,renorm_GF,ans2
max_diis_GW = 1
lin_GW = .false.
eta_GW = 0d0
reg_GW = .false.
TDA_W = .false.
read(1,*) if(ans1 == 'T') lin_GF = .true.
read(1,*) maxSCF_GW,thresh_GW,max_diis_GW,ans1,eta_GW,ans2,ans3 if(ans2 == 'T') reg_GF = .true.
if(ans1 == 'T') lin_GW = .true. ! Read GW options
if(ans2 == 'T') TDA_W = .true.
if(ans3 == 'T') reg_GW = .true.
! Read GT options maxSCF_GW = 64
thresh_GW = 1d-5
max_diis_GW = 1
lin_GW = .false.
eta_GW = 0d0
reg_GW = .false.
TDA_W = .false.
maxSCF_GT = 64 read(1,*)
thresh_GT = 1d-5 read(1,*) maxSCF_GW,thresh_GW,max_diis_GW,ans1,eta_GW,ans2,ans3
max_diis_GT = 1
lin_GT = .false.
eta_GT = 0d0
reg_GT = .false.
TDA_T = .false.
read(1,*) if(ans1 == 'T') lin_GW = .true.
read(1,*) maxSCF_GT,thresh_GT,max_diis_GT,ans1,eta_GT,ans2,ans3 if(ans2 == 'T') TDA_W = .true.
if(ans3 == 'T') reg_GW = .true.
if(ans1 == 'T') lin_GT = .true. ! Read GT options
if(ans2 == 'T') TDA_T = .true.
if(ans3 == 'T') reg_GT = .true.
! Options for adiabatic connection maxSCF_GT = 64
thresh_GT = 1d-5
max_diis_GT = 1
lin_GT = .false.
eta_GT = 0d0
reg_GT = .false.
TDA_T = .false.
doACFDT = .false. read(1,*)
exchange_kernel = .false. read(1,*) maxSCF_GT,thresh_GT,max_diis_GT,ans1,eta_GT,ans2,ans3
doXBS = .false.
read(1,*) if(ans1 == 'T') lin_GT = .true.
read(1,*) ans1,ans2,ans3 if(ans2 == 'T') TDA_T = .true.
if(ans3 == 'T') reg_GT = .true.
if(ans1 == 'T') doACFDT = .true. ! Options for adiabatic connection
if(ans2 == 'T') exchange_kernel = .true.
if(ans3 == 'T') doXBS = .true.
! Options for dynamical BSE calculations doACFDT = .false.
exchange_kernel = .false.
doXBS = .false.
dophBSE = .false. read(1,*)
dophBSE2 = .false. read(1,*) ans1,ans2,ans3
doppBSE = .false.
dBSE = .false.
dTDA = .true.
read(1,*) if(ans1 == 'T') doACFDT = .true.
read(1,*) ans1,ans2,ans3,ans4,ans5 if(ans2 == 'T') exchange_kernel = .true.
if(ans3 == 'T') doXBS = .true.
if(ans1 == 'T') dophBSE = .true. ! Options for dynamical BSE calculations
if(ans2 == 'T') dophBSE2 = .true.
if(ans3 == 'T') doppBSE = .true.
if(ans4 == 'T') dBSE = .true.
if(ans5 == 'F') dTDA = .false.
! Close file with options dophBSE = .false.
dophBSE2 = .false.
doppBSE = .false.
dBSE = .false.
dTDA = .true.
read(1,*)
read(1,*) ans1,ans2,ans3,ans4,ans5
if(ans1 == 'T') dophBSE = .true.
if(ans2 == 'T') dophBSE2 = .true.
if(ans3 == 'T') doppBSE = .true.
if(ans4 == 'T') dBSE = .true.
if(ans5 == 'F') dTDA = .false.
endif
! Close file with options
close(unit=1) close(unit=1)
end subroutine end subroutine

90
src/make_ninja.py
View File

@ -1,6 +1,9 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
import os import os
import sys import sys
import subprocess
DEBUG=False DEBUG=False
try: try:
@ -20,41 +23,28 @@ if "QUACK_ROOT" not in os.environ:
QUACK_ROOT=os.environ["QUACK_ROOT"] QUACK_ROOT=os.environ["QUACK_ROOT"]
if not DEBUG:
compile_gfortran_mac = """ def check_compiler_exists(compiler):
"""Check if a compiler exists on the system."""
try:
# Try to run the compiler with the --version flag to check its existence
subprocess.run([compiler, '--version'], check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
return True
except (subprocess.CalledProcessError, FileNotFoundError):
return False
compile_gfortran_mac = """
FC = gfortran FC = gfortran
AR = libtool -static -o AR = libtool -static -o
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -g -Wall -Wno-unused-variable -Wno-unused -Wno-unused-dummy-argument -O3 FFLAGS = -I$IDIR -J$IDIR -fbacktrace -g -Wall -Wno-unused-variable -Wno-unused -Wno-unused-dummy-argument -Wuninitialized -Wmaybe-uninitialized -O3 -march=native
CC = gcc CC = gcc
CXX = g++ CXX = g++
LAPACK=-lblas -llapack LAPACK=-lblas -llapack
STDCXX=-lc++ STDCXX=-lc++
FIX_ORDER_OF_LIBS=
"""
compile_gfortran_linux = """
FC = gfortran
AR = ar crs
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -g -Wall -Wno-unused -Wno-unused-dummy-argument -O3
CC = gcc
CXX = g++
LAPACK=-lblas -llapack
STDCXX=-lstdc++
FIX_ORDER_OF_LIBS=-Wl,--start-group FIX_ORDER_OF_LIBS=-Wl,--start-group
""" """
compile_ifort_linux = """ compile_gfortran_mac_debug = """
FC = ifort -mkl=parallel -qopenmp
AR = ar crs
FFLAGS = -I$IDIR -g -Ofast -traceback
CC = icc
CXX = icpc
LAPACK=
STDCXX=-lstdc++
FIX_ORDER_OF_LIBS=-Wl,--start-group
"""
else:
compile_gfortran_mac = """
FC = gfortran FC = gfortran
AR = libtool -static -o AR = libtool -static -o
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -Wno-unused-variable -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -Wno-unused-variable -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
@ -65,7 +55,7 @@ STDCXX=-lc++
FIX_ORDER_OF_LIBS= FIX_ORDER_OF_LIBS=
""" """
compile_gfortran_linux = """ compile_gfortran_linux_debug = """
FC = gfortran FC = gfortran
AR = ar crs AR = ar crs
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
@ -76,26 +66,50 @@ STDCXX=-lstdc++
FIX_ORDER_OF_LIBS=-Wl,--start-group FIX_ORDER_OF_LIBS=-Wl,--start-group
""" """
compile_olympe = """
FC = ifort -mkl=parallel -qopenmp
if sys.platform == "Darwin":
if DEBUG:
compiler = compile_gfortran_mac_debug
else:
compiler = compile_gfortran_mac
elif sys.platform == "Linux" or os.path.exists('/proc/version'):
if DEBUG:
compiler = compile_gfortran_linux_debug
else:
if check_compiler_exists('ifort'):
compiler = """
FC = ifort -qmkl=parallel -qopenmp
AR = ar crs AR = ar crs
FFLAGS = -I$IDIR -Ofast -traceback -xCORE-AVX512 FFLAGS = -I$IDIR -module $IDIR -traceback -g -Ofast -xHost
CC = icc CC = icc
CXX = icpc CXX = icpc
LAPACK= LAPACK=
STDCXX=-lstdc++ STDCXX=-lstdc++
FIX_ORDER_OF_LIBS=-Wl,--start-group FIX_ORDER_OF_LIBS=-Wl,--start-group
""" """
elif check_compiler_exists('gfortran'):
compiler = """
FC = gfortran -fopenmp
AR = ar crs
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -g -Wall -Wno-unused-variable -Wno-unused -Wno-unused-dummy-argument -Wuninitialized -Wmaybe-uninitialized -O3 -march=native
CC = gcc
CXX = g++
LAPACK=-lblas -llapack
STDCXX=-lstdc++
FIX_ORDER_OF_LIBS=-Wl,--start-group
"""
else:
raise RuntimeError("Neither ifort nor gfortran compilers were found on this system.")
if sys.platform in ["linux", "linux2"]:
# compiler = compile_gfortran_linux
compiler = compile_ifort_linux
# compiler = compile_olympe
elif sys.platform == "darwin":
compiler = compile_gfortran_mac
else: else:
print("Unknown platform. Only Linux and Darwin are supported.")
sys.exit(-1) print("Unknown platform. Only Linux and Darwin are supported.")
sys.exit(-1)
header = """# header = """#
# This file was automatically generated. Do not modify this file. # This file was automatically generated. Do not modify this file.

32
src/utils/read_basis_pyscf.f90
View File

@ -1,4 +1,4 @@
subroutine read_basis_pyscf(nBas,nO,nV) subroutine read_basis_pyscf(working_dir,nBas,nO,nV)
! Read basis set information from PySCF ! Read basis set information from PySCF
@ -7,21 +7,37 @@ subroutine read_basis_pyscf(nBas,nO,nV)
! Input variables ! Input variables
integer,intent(in) :: nO(nspin) integer,intent(in) :: nO(nspin)
character(len=256),intent(in) :: working_dir
! Local variables
! Output variables ! Output variables
integer,intent(out) :: nV(nspin) integer,intent(out) :: nV(nspin)
integer,intent(out) :: nBas integer,intent(out) :: nBas
! Local variables
integer :: status
character(len=256) :: file_path
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Primary basis set information ! Primary basis set information
!------------------------------------------------------------------------ !------------------------------------------------------------------------
open(unit=3,file='int/nBas.dat') file_path = trim(working_dir) // '/int/nBas.dat'
read(3, *) nBas open(unit=3, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
read(3, *) nBas
endif
close(unit=3) close(unit=3)
! write(*,'(A38)') '--------------------------------------' ! write(*,'(A38)') '--------------------------------------'

98
src/utils/read_dipole_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine read_dipole_integrals(nBas,R) subroutine read_dipole_integrals(working_dir,nBas,R)
! Read one-electron integrals related to dipole moment from files ! Read one-electron integrals related to dipole moment from files
@ -8,6 +8,7 @@ subroutine read_dipole_integrals(nBas,R)
! Input variables ! Input variables
integer,intent(in) :: nBas integer,intent(in) :: nBas
character(len=256),intent(in) :: working_dir
! Local variables ! Local variables
@ -19,36 +20,83 @@ subroutine read_dipole_integrals(nBas,R)
double precision,intent(out) :: R(nBas,nBas,ncart) double precision,intent(out) :: R(nBas,nBas,ncart)
integer :: status, ios
character(len=256) :: file_path
! Open file with integrals ! Open file with integrals
open(unit=21,file='int/x.dat')
open(unit=22,file='int/y.dat')
open(unit=23,file='int/z.dat')
! Read (x,y,z) integrals
R(:,:,:) = 0d0 R(:,:,:) = 0d0
do file_path = trim(working_dir) // '/int/x.dat'
read(21,*,end=21) mu,nu,Dip open(unit=21, file=file_path, status='old', action='read', iostat=status)
R(mu,nu,1) = Dip
R(nu,mu,1) = Dip
end do
21 close(unit=21)
do if(status /= 0) then
read(22,*,end=22) mu,nu,Dip
R(mu,nu,2) = Dip print *, "Error opening file: ", file_path
R(nu,mu,2) = Dip stop
end do
22 close(unit=22) else
do
read(21,*,iostat=ios) mu,nu,Dip
if(ios /= 0) exit
R(mu,nu,1) = Dip
R(nu,mu,1) = Dip
end do
endif
close(unit=21)
! ---
file_path = trim(working_dir) // '/int/y.dat'
open(unit=22, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
do
read(22,*,iostat=ios) mu,nu,Dip
if(ios /= 0) exit
R(mu,nu,2) = Dip
R(nu,mu,2) = Dip
end do
endif
close(unit=22)
! ---
file_path = trim(working_dir) // '/int/z.dat'
open(unit=23, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
do
read(23,*,iostat=ios) mu,nu,Dip
if(ios /= 0) exit
R(mu,nu,3) = Dip
R(nu,mu,3) = Dip
end do
endif
close(unit=23)
! ---
do
read(23,*,end=23) mu,nu,Dip
R(mu,nu,3) = Dip
R(nu,mu,3) = Dip
end do
23 close(unit=23)
! Print results ! Print results
if(debug) then if(debug) then

85
src/utils/read_geometry.f90
View File

@ -1,10 +1,14 @@
subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc) subroutine read_geometry(working_dir,nNuc,ZNuc,rNuc,ENuc)
! Read molecular geometry ! Read molecular geometry
implicit none implicit none
include 'parameters.h' include 'parameters.h'
! Input variables
character(len=256),intent(in) :: working_dir
! Ouput variables ! Ouput variables
integer,intent(in) :: nNuc integer,intent(in) :: nNuc
@ -20,43 +24,70 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
double precision,intent(out) :: ZNuc(nNuc),rNuc(nNuc,ncart),ENuc double precision,intent(out) :: ZNuc(nNuc),rNuc(nNuc,ncart),ENuc
integer :: status
character(len=256) :: file_path
! Open file with geometry specification ! Open file with geometry specification
open(unit=10,file='input/molecule') file_path = trim(working_dir) // '/input/molecule'
open(unit=11,file='input/molecule.xyz') open(unit=10, file=file_path, status='old', action='read', iostat=status)
! Read geometry and create xyz file for integrals if(status /= 0) then
read(10,*) print *, "Error opening file: ", file_path
read(10,*) stop
read(10,*)
write(11,'(I3)') nNuc else
write(11,*)
do i=1,nNuc ! Read geometry and create xyz file for integrals
read(10,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3) open(unit=11,file=trim(working_dir) // '/input/molecule.xyz')
write(11,'(A3,1X,3F16.10)') El,rNuc(i,1)*BoToAn,rNuc(i,2)*BoToAn,rNuc(i,3)*BoToAn
ZNuc(i) = dble(element_number(El))
end do
! Compute nuclear repulsion energy read(10,*)
read(10,*)
read(10,*)
write(11,'(I3)') nNuc
write(11,*)
do i=1,nNuc
read(10,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3)
write(11,'(A3,1X,3F16.10)') El,rNuc(i,1)*BoToAn,rNuc(i,2)*BoToAn,rNuc(i,3)*BoToAn
ZNuc(i) = dble(element_number(El))
end do
close(unit=11)
endif
close(unit=10)
! ---
file_path = trim(working_dir) // '/int/ENuc.dat'
open(unit=3, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
read(3,*) ENuc
endif
ENuc = 0
open(unit=3,file='int/ENuc.dat')
read(3,*) ENuc
close(unit=3) close(unit=3)
! do i=1,nNuc-1 ! Compute nuclear repulsion energy
! do j=i+1,nNuc !ENuc = 0.d0
! RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2 !do i=1,nNuc-1
! ENuc = ENuc + ZNuc(i)*ZNuc(j)/(AntoBo*sqrt(RAB)) ! do j=i+1,nNuc
! end do ! RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2
! end do ! ENuc = ENuc + ZNuc(i)*ZNuc(j)/(AntoBo*sqrt(RAB))
! end do
!end do
! Close file with geometry specification
close(unit=10)
close(unit=11)
! Print geometry ! Print geometry
write(*,'(A28)') '------------------' write(*,'(A28)') '------------------'

104
src/utils/read_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine read_integrals(nBas_AOs, S, T, V, Hc, G) subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
! Read one- and two-electron integrals from files ! Read one- and two-electron integrals from files
@ -8,6 +8,7 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
! Input variables ! Input variables
integer,intent(in) :: nBas_AOs integer,intent(in) :: nBas_AOs
character(len=256),intent(in) :: working_dir
! Local variables ! Local variables
@ -24,6 +25,9 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
double precision,intent(out) :: Hc(nBas_AOs,nBas_AOs) double precision,intent(out) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
integer :: status, ios
character(len=256) :: file_path
! Open file with integrals ! Open file with integrals
debug = .false. debug = .false.
@ -32,46 +36,67 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
print*, 'Scaling integrals by ',lambda print*, 'Scaling integrals by ',lambda
open(unit=8 ,file='int/Ov.dat')
open(unit=9 ,file='int/Kin.dat')
open(unit=10,file='int/Nuc.dat')
open(unit=21,file='int/x.dat') ! ---
open(unit=22,file='int/y.dat')
open(unit=23,file='int/z.dat')
! Read overlap integrals ! Read overlap integrals
file_path = trim(working_dir) // '/int/Ov.dat'
open(unit=8, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
S(:,:) = 0d0
do
read(8,*,iostat=ios) mu,nu,Ov
if(ios /= 0) exit
S(mu,nu) = Ov
S(nu,mu) = Ov
end do
endif
close(unit=8)
S(:,:) = 0d0 ! ---
do
read(8,*,end=8) mu,nu,Ov
S(mu,nu) = Ov
S(nu,mu) = Ov
end do
8 close(unit=8)
! Read kinetic integrals ! Read kinetic integrals
file_path = trim(working_dir) // '/int/Kin.dat'
open(unit=9, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
T(:,:) = 0d0
do
read(9,*,iostat=ios) mu,nu,Kin
if(ios /= 0) exit
T(mu,nu) = Kin
T(nu,mu) = Kin
end do
endif
close(unit=9)
T(:,:) = 0d0 ! ---
do
read(9,*,end=9) mu,nu,Kin
T(mu,nu) = Kin
T(nu,mu) = Kin
end do
9 close(unit=9)
! Read nuclear integrals ! Read nuclear integrals
file_path = trim(working_dir) // '/int/Nuc.dat'
open(unit=10, file=file_path, status='old', action='read', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
V(:,:) = 0d0
do
read(10,*,iostat=ios) mu,nu,Nuc
if(ios /= 0) exit
V(mu,nu) = Nuc
V(nu,mu) = Nuc
end do
endif
close(unit=10)
V(:,:) = 0d0 ! ---
do
read(10,*,end=10) mu,nu,Nuc
V(mu,nu) = Nuc
V(nu,mu) = Nuc
end do
10 close(unit=10)
! Define core Hamiltonian
! Define core Hamiltonian
Hc(:,:) = T(:,:) + V(:,:) Hc(:,:) = T(:,:) + V(:,:)
! Read 2e-integrals ! Read 2e-integrals
@ -94,9 +119,16 @@ subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
! 11 close(unit=11) ! 11 close(unit=11)
! binary file ! binary file
open(unit=11, file='int/ERI.bin', form='unformatted', access='stream') file_path = trim(working_dir) // '/int/ERI.bin'
read(11) G open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=status)
close(11) if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
read(11) G
endif
close(unit=11)
! Print results ! Print results

72
src/utils/read_molecule.f90
View File

@ -1,4 +1,4 @@
subroutine read_molecule(nNuc,nO,nC,nR) subroutine read_molecule(working_dir,nNuc,nO,nC,nR)
! Read number of atoms and number of electrons ! Read number of atoms and number of electrons
@ -11,6 +11,10 @@ subroutine read_molecule(nNuc,nO,nC,nR)
integer :: nCore integer :: nCore
integer :: nRyd integer :: nRyd
! Input variables
character(len=256),intent(in) :: working_dir
! Output variables ! Output variables
integer,intent(out) :: nNuc integer,intent(out) :: nNuc
@ -18,45 +22,57 @@ subroutine read_molecule(nNuc,nO,nC,nR)
integer,intent(out) :: nC(nspin) integer,intent(out) :: nC(nspin)
integer,intent(out) :: nR(nspin) integer,intent(out) :: nR(nspin)
integer :: status
character(len=256) :: file_path
! Open file with geometry specification ! Open file with geometry specification
open(unit=1,file='input/molecule') file_path = trim(working_dir) // '/input/molecule'
open(unit=1, file=file_path, status='old', action='read', iostat=status)
! Read number of atoms and number of electrons if(status /= 0) then
read(1,*) print *, "Error opening file: ", file_path
read(1,*) nNuc,nO(1),nO(2),nCore,nRyd stop
if(mod(nCore,2) /= 0 .or. mod(nRyd,2) /= 0) then else
print*, 'The number of core and Rydberg electrons must be even!' ! Read number of atoms and number of electrons
stop
end if read(1,*)
read(1,*) nNuc,nO(1),nO(2),nCore,nRyd
nC(:) = nCore/2 if(mod(nCore,2) /= 0 .or. mod(nRyd,2) /= 0) then
nR(:) = nRyd/2
! Print results print*, 'The number of core and Rydberg electrons must be even!'
stop
write(*,'(A28)') '----------------------' end if
write(*,'(A28,1X,I16)') 'Number of atoms',nNuc
write(*,'(A28)') '----------------------'
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nO(1)
write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nO(2)
write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nO(:))
write(*,'(A28)') '----------------------'
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of core electrons',sum(nC(:))
write(*,'(A28,1X,I16)') 'Number of Rydberg electrons',sum(nR(:))
write(*,'(A28)') '----------------------'
write(*,*)
! Close file with geometry specification nC(:) = nCore/2
nR(:) = nRyd/2
! Print results
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of atoms',nNuc
write(*,'(A28)') '----------------------'
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nO(1)
write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nO(2)
write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nO(:))
write(*,'(A28)') '----------------------'
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of core electrons',sum(nC(:))
write(*,'(A28,1X,I16)') 'Number of Rydberg electrons',sum(nR(:))
write(*,'(A28)') '----------------------'
write(*,*)
endif
! Close file with geometry specification
close(unit=1) close(unit=1)
end subroutine end subroutine