LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

Merge pull request #24 from LCPQ/master

Browse Source 
Merge with eginer
This commit is contained in:
Thomas Applencourt 2015-07-03 12:24:01 +02:00
commit e91c676dd6
12 changed files with 617 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
plugins/MRCC/davidson.irp.f
View File

@ -389,16 +389,17 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
Vt = 0.d0
!$OMP DO SCHEDULE(guided)
do i=1,n
idx(0) = i
call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
do jj=1,idx(0)
j = idx(jj)
if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
! idx(0) = i
! call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
! do jj=1,idx(0)
! j = idx(jj)
! if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
do j = 1, i-1
call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
hij = hij + delta_ij(j,i,istate)
vt (i) = vt (i) + hij*u_0(j)
vt (j) = vt (j) + hij*u_0(i)
endif
! endif
enddo
enddo
!$OMP END DO

2
plugins/MRCC/mrcc.irp.f
View File

@ -48,7 +48,7 @@ subroutine run_mrcc
E_new = 0.d0
delta_E = 1.d0
iteration = 0
do while (delta_E > 1.d-8)
do while (delta_E > 1.d-10)
iteration += 1
print *, '==========================='
print *, 'MRCC Iteration', iteration

4
plugins/MRCC/mrcc_dress.irp.f
View File

@ -138,7 +138,11 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro
do i_state=1,N_states
delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),i_state) += dIa_hla(i_state,k_sd)
delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),i_state) += dIa_hla(i_state,k_sd)
if(dabs(psi_cas_coef(i_I,i_state)).ge.5.d-5)then
delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state)
else
delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) = 0.d0
endif
enddo
enddo
call omp_unset_lock( psi_cas_lock(i_I) )

37
plugins/MRCC/mrcc_utils.irp.f
View File

@ -1,35 +1,28 @@
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states,psi_det_size) ]
implicit none
BEGIN_DOC
! cm/<Psi_0|H|D_m>
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
END_DOC
integer :: i,k
double precision :: ihpsi(N_states), hij(N_states)
double precision :: ihpsi(N_states), hii
do i=1,N_det_non_cas
call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coef, N_int, N_det_cas, &
size(psi_cas_coef,1), n_states, ihpsi)
call i_h_j(psi_non_cas(1,1,i),psi_non_cas(1,1,i),N_int,hij)
call i_h_j(psi_non_cas(1,1,i),psi_non_cas(1,1,i),N_int,hii)
do k=1,N_states
lambda_pert(k,i) = 1d0 / (CI_electronic_energy(k)-hij(k))
lambda_pert(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii)
lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k)
if ((lambda_mrcc(k,i)/lambda_pert(k,i))<0.d0 .or. (lambda_mrcc(k,i)/lambda_pert(k,i))>4.d0) then
lambda_mrcc(k,i) = lambda_pert(k,i)
else
if ((lambda_mrcc(k,i)/lambda_pert(k,i))<0.1d0 .or. (lambda_mrcc(k,i)/lambda_pert(k,i))>=0d0) then
lambda_mrcc(k,i) = lambda_mrcc(k,i)*((cos((lambda_mrcc(k,i)/lambda_pert(k,i))*3.141592653589793d0/0.1d0+3.141592653589793d0)+1d0)/2.d0) &
+ lambda_pert(k,i)*(1.d0-((cos((lambda_mrcc(k,i)/lambda_pert(k,i))*3.141592653589793d0/0.1d0+3.141592653589793d0)+1.d0)/2.d0))
elseif ((lambda_mrcc(k,i)/lambda_pert(k,i))<=4.0d0 .or. (lambda_mrcc(k,i)/lambda_pert(k,i))>2.0d0) then
lambda_mrcc(k,i) = lambda_mrcc(k,i)*(1.d0-(cos(abs(2.d0-(lambda_mrcc(k,i)/lambda_pert(k,i)))*3.141592653589793d0/2.0d0+3.141592653589793d0)+1.d0)/2d0) &
+ lambda_pert(k,i)*((cos(abs(2.d0-(lambda_mrcc(k,i)/lambda_pert(k,i)))*3.141592653589793d0/2.0d0+3.141592653589793d0)+1.d0)/2.d0)
else
lambda_mrcc(k,i) = lambda_mrcc(k,i)
endif
if (dabs(ihpsi(k)).le.1.d-3) then
lambda_mrcc(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii)
icount_manu = icount_manu+1
cycle
endif
enddo
enddo
END_PROVIDER
@ -71,6 +64,16 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
enddo
enddo
endif
do i = 1, N_det
do j = 1, N_det
do m = 1, N_states
if(isnan(delta_ij(j,i,m)))then
delta_ij(j,i,m) = 0.d0
endif
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]

1
plugins/loc_cele/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
AO_Basis Electrons Ezfio_files MO_Basis Nuclei Utils

28
plugins/loc_cele/README.rst Normal file
View File

@ -0,0 +1,28 @@
===============
loc_cele Module
===============
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. by the `update_README.py` script.
`loc_rasorb <http://github.com/LCPQ/quantum_package/tree/master/src/loc_cele/loc_cele.irp.f#L1>`_
Undocumented
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. by the `update_README.py` script.
.. image:: tree_dependency.png
* `AO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/AO_Basis>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_

163
plugins/loc_cele/loc.f Normal file
View File

@ -0,0 +1,163 @@
c************************************************************************
subroutine maxovl(n,m,s,t,w)
C
C This subprogram contains an iterative procedure to find the
C unitary transformation of a set of n vectors which maximizes
C the sum of their square overlaps with a set of m reference
C vectors (m.le.n)
C
C S: overlap matrix <ref|vec>
C T: rotation matrix
C W: new overlap matrix
C
C
implicit real*8(a-h,o-y),logical*1(z)
parameter (id1=300)
dimension s(id1,id1),t(id1,id1),w(id1,id1)
data small/1.d-6/
zprt=.true.
niter=100
conv=1.d-8
write (6,5) n,m,conv
5 format (//5x,'Unitary transformation of',i3,' vectors'/
* 5x,'following the principle of maximum overlap with a set of',
* i3,' reference vectors'/5x,'required convergence on rotation ',
* 'angle =',f13.10///5x,'Starting overlap matrix'/)
do 6 i=1,m
write (6,145) i
6 write (6,150) (s(i,j),j=1,n)
8 mm=m-1
if (m.lt.n) mm=m
iter=0
do 20 j=1,n
do 16 i=1,n
t(i,j)=0.d0
16 continue
do 18 i=1,m
18 w(i,j)=s(i,j)
20 t(j,j)=1.d0
sum=0.d0
do 10 i=1,m
sum=sum+s(i,i)*s(i,i)
10 continue
sum=sum/m
if (zprt) write (6,12) sum
12 format (//5x,'Average square overlap =',f10.6)
if (n.eq.1) goto 100
last=n
j=1
21 if (j.ge.last) goto 30
sum=0.d0
do 22 i=1,n
22 sum=sum+s(i,j)*s(i,j)
if (sum.gt.small) goto 28
do 24 i=1,n
sij=s(i,j)
s(i,j)=-s(i,last)
s(i,last)=sij
tij=t(i,j)
t(i,j)=-t(i,last)
t(i,last)=tij
24 continue
last=last-1
goto 21
28 j=j+1
goto 21
30 iter=iter+1
imax=0
jmax=0
dmax=0.d0
amax=0.d0
do 60 i=1,mm
ip=i+1
do 50 j=ip,n
a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
b=-s(i,i)*s(i,j)
if (j.gt.m) goto 31
a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
b=b+s(j,i)*s(j,j)
31 b=b+b
if (a.eq.0.d0) goto 32
ba=b/a
if (dabs(ba).gt.small) goto 32
if (a.gt.0.d0) goto 33
tang=-0.5d0*ba
cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine
goto 34
32 tang=0.d0
if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine
goto 34
33 cosine=0.d0
sine=1.d0
34 delta=sine*(a*sine+b*cosine)
if (zprt.and.delta.lt.0.d0) write (6,71) i,j,a,b,sine,cosine,delta
do 35 k=1,m
p=s(k,i)*cosine-s(k,j)*sine
q=s(k,i)*sine+s(k,j)*cosine
s(k,i)=p
35 s(k,j)=q
do 40 k=1,n
p=t(k,i)*cosine-t(k,j)*sine
q=t(k,i)*sine+t(k,j)*cosine
t(k,i)=p
t(k,j)=q
40 continue
45 d=dabs(sine)
if (d.le.amax) goto 50
imax=i
jmax=j
amax=d
dmax=delta
50 continue
60 continue
if (zprt) write (6,70) iter,amax,imax,jmax,dmax
70 format (' iter=',i4,' largest rotation=',f12.8,
* ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
if (amax.lt.conv) goto 100
if (iter.lt.niter) goto 30
write (6,80)
write (6,*) 'niter=',niter
80 format (//5x,'*** maximum number of cycles exceeded ',
* 'in subroutine maxovl ***'//)
stop
100 continue
do 120 j=1,n
if (s(j,j).gt.0.d0) goto 120
do 105 i=1,m
105 s(i,j)=-s(i,j)
do 110 i=1,n
110 t(i,j)=-t(i,j)
120 continue
sum=0.d0
do 125 i=1,m
125 sum=sum+s(i,i)*s(i,i)
sum=sum/m
do 122 i=1,m
do 122 j=1,n
sw=s(i,j)
s(i,j)=w(i,j)
122 w(i,j)=sw
if (.not.zprt) return
write (6,12) sum
write (6,130)
130 format (//5x,'transformation matrix')
do 140 i=1,n
write (6,145) i
140 write (6,150) (t(i,j),j=1,n)
145 format (i8)
150 format (2x,10f12.8)
write (6,160)
160 format (//5x,'new overlap matrix'/)
do 170 i=1,m
write (6,145) i
170 write (6,150) (w(i,j),j=1,n)
return
end

333
plugins/loc_cele/loc_cele.irp.f Normal file
View File

@ -0,0 +1,333 @@
program loc_rasorb
implicit none
!
! This program performs a localization of the active orbitals
! of a CASSCF wavefunction, reading the orbitals from a RASORB
! file of molcas.
!
! id1=max number of MO in a given symmetry.
!
integer id1
parameter (id1=300)
character*1 jobz,uplo
character*64 file1,file2
character*72 string(id1,8),cdum
double precision :: cmo(id1,id1,1),cmoref(id1,id1,1),newcmo(id1,id1,1)
double precision ::s(id1,id1,1),dum,ddum(id1,id1),ovl(id1,id1)
double precision :: w(id1),work(3*id1),t(id1,id1),wi(id1,id1)
integer n,i,j,k,l,nmo(8),isym,nsym,idum,nrot(8),irot(id1,8)
integer ipiv(id1),info,lwork
logical *1 z54
print*,'passed the first copy'
z54=.false.
!Read the name of the RasOrb file
print*,'Entering in the loc program'
! read(5,*) z54
print*,'before = '
accu_norm = 0.d0
do i =1,mo_tot_num
accu_norm += dabs(mo_overlap(i,i))
enddo
print*,'accu_norm = ',accu_norm
nsym = 1
nmo(1) = mo_tot_num
print*,'nmo(1) = ',nmo(1)
cmo = 0.d0
do isym=1,nsym
do i=1,nmo(isym)
do j = 1, ao_num
cmo(j,i,isym) = mo_coef(j,i)
enddo
enddo
enddo
print*,'passed the first copy'
do isym=1,nsym
do j=1,mo_tot_num
do i=1,ao_num
newcmo(i,j,isym)=cmo(i,j,isym)
enddo
enddo
enddo
print*,'passed the copy'
nrot(1) = 6 ! number of orbitals to be localized
integer :: index_rot(1000,1)
cmoref = 0.d0
! Definition of the index of the MO to be rotated
irot(1,1) = 20 ! the first mo to be rotated is the 19 th MO
irot(2,1) = 21 ! the first mo to be rotated is the 20 th MO
irot(3,1) = 22 ! etc....
irot(4,1) = 23 !
irot(5,1) = 24 !
irot(6,1) = 25 !
! you define the guess vectors that you want
! the new MO to be close to
! cmore(i,j,1) = < AO_i | guess_vector_MO(j) >
! i goes from 1 to ao_num
! j goes from 1 to nrot(1)
! Here you must go to the GAMESS output file
! where the AOs are listed and explicited
! From the basis of this knowledge you can build your
! own guess vectors for the MOs
! The new MOs are provided in output
! in the same order than the guess MOs
cmoref(3,1,1) = 1.d0 !
cmoref(12,1,1) = 1.d0 !
cmoref(21,2,1) = 1.d0 !
cmoref(30,2,1) = 1.d0 !
cmoref(39,3,1) = 1.d0 !
cmoref(48,3,1) = 1.d0 !
cmoref(3,4,1) = 1.d0 !
cmoref(12,4,1) =-1.d0 !
cmoref(21,5,1) = 1.d0 !
cmoref(30,5,1) =-1.d0 !
cmoref(39,6,1) = 1.d0 !
cmoref(48,6,1) =-1.d0 !
print*,'passed the definition of the referent vectors '
!Building the S (overlap) matrix in the AO basis.
do isym=1,nsym
if (nrot(isym).eq.0) cycle
do i=1,ao_num
s(i,i,isym)=1.d0
do j=1,ao_num
if (i.ne.j) s(i,j,isym)=0.d0
ddum(i,j)=0.d0
do k=1,nmo(isym)
ddum(i,j)=ddum(i,j)+cmo(i,k,isym)*cmo(j,k,isym)
enddo
enddo
enddo
call dgesv(ao_num,ao_num,ddum,id1,ipiv,s(1,1,isym),id1,info)
if (info.ne.0) then
write (6,*) 'Something wrong in dgsev',isym
stop
endif
enddo
!Now big loop over symmetry
do isym=1,nsym
if (nrot(isym).eq.0) cycle
write (6,*)
write (6,*)
write (6,*)
write (6,*) 'WORKING ON SYMMETRY',isym
write (6,*)
!Compute the overlap matrix <ref|vec>
! do i=1,nmo(isym)
do i=1,ao_num
do j=1,nrot(isym)
ddum(i,j)=0.d0
do k=1,ao_num
ddum(i,j)=ddum(i,j)+s(i,k,isym)*cmo(k,irot(j,isym),isym)
enddo
enddo
enddo
do i=1,nrot(isym)
do j=1,nrot(isym)
ovl(i,j)=0.d0
do k=1,ao_num
! do k=1,mo_tot_num
ovl(i,j)=ovl(i,j)+cmoref(k,i,isym)*ddum(k,j)
enddo
enddo
enddo
call maxovl(nrot(isym),nrot(isym),ovl,t,wi)
do i=1,nrot(isym)
do j=1,ao_num
write (6,*) 'isym,',isym,nrot(isym),nmo(isym)
newcmo(j,irot(i,isym),isym)=0.d0
do k=1,nrot(isym)
newcmo(j,irot(i,isym),isym)=newcmo(j,irot(i,isym),isym) + cmo(j,irot(k,isym),isym)*t(k,i)
enddo
enddo
enddo
! if(dabs(newcmo(3,19,1) - mo_coef(3,19)) .gt.1.d-10 )then
! print*,'Something wrong bitch !!'
! print*,'newcmo(3,19,1) = ',newcmo(3,19,1)
! print*,'mo_coef(3,19) = ',mo_coef(3,19)
! stop
! endif
enddo !big loop over symmetry
10 format (4E18.12)
! Now we copyt the newcmo into the mo_coef
mo_coef = 0.d0
do isym=1,nsym
do i=1,nmo(isym)
do j = 1, ao_num
mo_coef(j,i) = newcmo(j,i,isym)
enddo
enddo
enddo
! if(dabs(newcmo(3,19,1) - mo_coef(3,19)) .gt.1.d-10 )then
print*,'mo_coef(3,19)',mo_coef(3,19)
pause
! we say that it hase been touched, and valid and that everything that
! depends on mo_coef must not be reprovided
double precision :: accu_norm
touch mo_coef
print*,'after = '
accu_norm = 0.d0
do i =1,mo_tot_num
accu_norm += dabs(mo_overlap(i,i))
enddo
print*,'accu_norm = ',accu_norm
! We call the routine that saves mo_coef in the ezfio format
call save_mos
stop
end

BIN
plugins/loc_cele/tree_dependency.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 48 KiB

2
src/Determinants/davidson.irp.f
View File

@ -376,7 +376,7 @@ end
! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
END_DOC
davidson_criterion = 'residual'
davidson_threshold = 1.d-6
davidson_threshold = 1.d-9
END_PROVIDER
subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)

1
src/Determinants/occ_pattern.irp.f
View File

@ -256,6 +256,7 @@ subroutine make_s2_eigenfunction
integer :: N_det_new
integer, parameter :: bufsze = 1000
logical, external :: is_in_wavefunction
return
! !TODO DEBUG
! do i=1,N_det

54
src/Determinants/psi_cas.irp.f
View File

@ -13,6 +13,9 @@ use bitmasks
logical :: good
N_det_cas = 0
do i=1,N_det
do l = 1, N_states
psi_cas_coef(i,l) = 0.d0
enddo
do l=1,n_cas_bitmask
good = .True.
do k=1,N_int
@ -109,6 +112,57 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [double precision, H_matrix_cas, (N_det_cas,N_det_cas)]
implicit none
integer :: i,j
double precision :: hij
do i = 1, N_det_cas
do j = 1, N_det_cas
call i_H_j(psi_cas(1,1,i),psi_cas(1,1,j),N_int,hij)
H_matrix_cas(i,j) = hij
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, psi_coef_cas_diagonalized, (N_det_cas,N_states)]
&BEGIN_PROVIDER [double precision, psi_cas_energy_diagonalized, (N_states)]
implicit none
integer :: i,j
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
allocate (eigenvectors(size(H_matrix_cas,1),N_det_cas))
allocate (eigenvalues(N_det_cas))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_cas,size(H_matrix_cas,1),N_det_cas)
do i = 1, N_states
psi_cas_energy_diagonalized(i) = eigenvalues(i)
do j = 1, N_det_cas
psi_coef_cas_diagonalized(j,i) = eigenvectors(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, psi_cas_energy, (N_states)]
implicit none
integer :: i,j,k
double precision :: hij,norm,u_dot_v
psi_cas_energy = 0.d0
do k = 1, N_states
norm = 0.d0
do i = 1, N_det_cas
norm += psi_cas_coef(i,k) * psi_cas_coef(i,k)
do j = 1, N_det_cas
psi_cas_energy(k) += psi_cas_coef(i,k) * psi_cas_coef(j,k) * H_matrix_cas(i,j)
enddo
enddo
psi_cas_energy(k) = psi_cas_energy(k) /norm
enddo
END_PROVIDER