diff --git a/plugins/MRCC/davidson.irp.f b/plugins/MRCC/davidson.irp.f index 52f03ff3..57900082 100644 --- a/plugins/MRCC/davidson.irp.f +++ b/plugins/MRCC/davidson.irp.f @@ -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 diff --git a/plugins/MRCC/mrcc.irp.f b/plugins/MRCC/mrcc.irp.f index 30f6645b..8a35dde5 100644 --- a/plugins/MRCC/mrcc.irp.f +++ b/plugins/MRCC/mrcc.irp.f @@ -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 diff --git a/plugins/MRCC/mrcc_dress.irp.f b/plugins/MRCC/mrcc_dress.irp.f index 083cea5a..23c15d13 100644 --- a/plugins/MRCC/mrcc_dress.irp.f +++ b/plugins/MRCC/mrcc_dress.irp.f @@ -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) - 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) + 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) ) diff --git a/plugins/MRCC/mrcc_utils.irp.f b/plugins/MRCC/mrcc_utils.irp.f index 716d5ffc..b175d655 100644 --- a/plugins/MRCC/mrcc_utils.irp.f +++ b/plugins/MRCC/mrcc_utils.irp.f @@ -1,35 +1,28 @@ - -BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] + 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/ + ! cm/ 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 - 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) ] diff --git a/plugins/loc_cele/NEEDED_CHILDREN_MODULES b/plugins/loc_cele/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..a23aa5ae --- /dev/null +++ b/plugins/loc_cele/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +AO_Basis Electrons Ezfio_files MO_Basis Nuclei Utils diff --git a/plugins/loc_cele/README.rst b/plugins/loc_cele/README.rst new file mode 100644 index 00000000..431dae8a --- /dev/null +++ b/plugins/loc_cele/README.rst @@ -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 `_ + 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 `_ +* `Electrons `_ +* `Ezfio_files `_ +* `MO_Basis `_ +* `Nuclei `_ +* `Utils `_ + diff --git a/plugins/loc_cele/loc.f b/plugins/loc_cele/loc.f new file mode 100644 index 00000000..e2439b7f --- /dev/null +++ b/plugins/loc_cele/loc.f @@ -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 +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 + diff --git a/plugins/loc_cele/loc_cele.irp.f b/plugins/loc_cele/loc_cele.irp.f new file mode 100644 index 00000000..5b9112ca --- /dev/null +++ b/plugins/loc_cele/loc_cele.irp.f @@ -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 + + + + + +! 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 diff --git a/plugins/loc_cele/tree_dependency.png b/plugins/loc_cele/tree_dependency.png new file mode 100644 index 00000000..244d9db0 Binary files /dev/null and b/plugins/loc_cele/tree_dependency.png differ diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index bdc979c4..e59fde4f 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -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) diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index a0fd4a3c..8027cf60 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/src/Determinants/occ_pattern.irp.f @@ -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 diff --git a/src/Determinants/psi_cas.irp.f b/src/Determinants/psi_cas.irp.f index 8ca081d6..304a2370 100644 --- a/src/Determinants/psi_cas.irp.f +++ b/src/Determinants/psi_cas.irp.f @@ -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 +