diff --git a/cdn/rot_den_mat.f b/cdn/rot_den_mat.f index faf5a8b7bfba1d96045095398b697bf08a7aeb7a..c5cd2eced45a3bd6ad62a481ac91f5c900ab9a97 100644 --- a/cdn/rot_den_mat.f +++ b/cdn/rot_den_mat.f @@ -15,6 +15,7 @@ c This subroutine rotates the direction of the magnetization of the c density matrix by multiplying with the unitary 2x2 spin rotation c matrix. --> U*rho*U^dagger c Philipp Kurz 2000-02-03 +c new method for improved stability (l_new=t) gb'19 c*********************************************************************** use m_constants @@ -28,11 +29,29 @@ C .. Scalar Arguments .. C .. C .. Local Scalars .. INTEGER ispin - REAL eps + REAL eps,r11n,r22n + COMPLEX r21n + LOGICAL l_new C .. C .. Local Arrays .. COMPLEX u2(2,2),rho(2,2),rhoh(2,2) C .. + l_new = .true. + + IF (l_new) THEN + + r11n = 0.5*(1.0+cos(beta))*rho11 - sin(beta)*real(rho21) + + + 0.5*(1.0-cos(beta))*rho22 + r22n = 0.5*(1.0-cos(beta))*rho11 + sin(beta)*real(rho21) + + + 0.5*(1.0+cos(beta))*rho22 + r21n = CMPLX(cos(alph),-sin(alph))*(sin(beta)*(rho11-rho22) + + + 2.0*(cos(beta)*real(rho21)-cmplx(0.0,aimag(rho21))))*0.5 + + rho11 = r11n + rho22 = r22n + rho21 = r21n + + ELSE eps = 1.0e-10 @@ -79,6 +98,8 @@ c---> are real. rho22 = real(rho(2,2)) rho21 = rho(2,1) + ENDIF + END SUBROUTINE rot_den_mat END MODULE m_rotdenmat diff --git a/eigen/hsvac.F90 b/eigen/hsvac.F90 index b1c3bbb53d3eb935322931cf14a147fff3272590..57dd3bd683f007b30472af2d07c14ee9a619fd12 100644 --- a/eigen/hsvac.F90 +++ b/eigen/hsvac.F90 @@ -29,7 +29,7 @@ CONTAINS TYPE(t_lapw),INTENT(IN) :: lapw TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_potden),INTENT(IN) :: v - CLASS(t_mat),INTENT(INOUT) :: hmat(:,:),smat(:,:) + CLASS(t_mat),INTENT(INOUT) :: hmat(:,:),smat(:,:) ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: jsp @@ -40,14 +40,14 @@ CONTAINS ! .. Local Scalars .. COMPLEX hij,sij,apw_lo,c_1 REAL d2,gz,sign,th,wronk - INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ipot,ii0,i0 + INTEGER i,i2,ii,jj,ik,j,jk,k,jspin,ii0,i0 INTEGER ivac,irec,imz,igvm2,igvm2i,s1,s2 INTEGER jspin1,jspin2,jmax,jsp_start,jsp_end INTEGER i_start,nc,nc_0 ! .. ! .. Local Arrays .. INTEGER:: nv2(input%jspins) - INTEGER kvac1(DIMENSION%nv2d,input%jspins),kvac2(DIMENSION%nv2d,input%jspins) + INTEGER kvac(2,DIMENSION%nv2d,input%jspins) INTEGER map2(DIMENSION%nvd,input%jspins) COMPLEX tddv(DIMENSION%nv2d,DIMENSION%nv2d),tduv(DIMENSION%nv2d,DIMENSION%nv2d) COMPLEX tudv(DIMENSION%nv2d,DIMENSION%nv2d),tuuv(DIMENSION%nv2d,DIMENSION%nv2d) @@ -66,16 +66,14 @@ CONTAINS nv2(jspin) = 0 k_loop:DO k = 1,lapw%nv(jspin) DO j = 1,nv2(jspin) - IF (lapw%k1(k,jspin).EQ.kvac1(j,jspin)& - .AND. lapw%k2(k,jspin).EQ.kvac2(j,jspin)) THEN + IF (all(lapw%gvec(1:2,k,jspin)==kvac(1:2,j,jspin))) THEN map2(k,jspin) = j CYCLE k_loop END IF ENDDO nv2(jspin) = nv2(jspin) + 1 IF (nv2(jspin)>DIMENSION%nv2d) CALL juDFT_error("hsvac:dimension%nv2d",calledby ="hsvac") - kvac1(nv2(jspin),jspin) = lapw%k1(k,jspin) - kvac2(nv2(jspin),jspin) = lapw%k2(k,jspin) + kvac(1:2,nv2(jspin),jspin) = lapw%gvec(1:2,k,jspin) map2(k,jspin) = nv2(jspin) ENDDO k_loop ENDDO @@ -86,20 +84,16 @@ CONTAINS s1=MIN(SIZE(hmat,1),jspin1) !in colinear case s1=1 DO jspin2=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin s2=MIN(SIZE(hmat,1),jspin2) !in colinear case s2=1 - ipot=3 - IF (jspin1==jspin2) ipot=jspin1 !---> get the wavefunctions and set up the tuuv, etc matrices - jspin=jsp - CALL vacfun(& - vacuum,DIMENSION,stars,& - jsp,input,noco,jspin1,jspin2,& - sym, cell,ivac,evac(1,1),lapw%bkpt,v%vacxy(:,:,ivac,ipot),v%vacz(:,:,:),kvac1,kvac2,nv2,& - tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk) + CALL vacfun(& + vacuum,stars,& + input,noco,jspin1,jspin2,& + sym, cell,ivac,evac,lapw%bkpt,v%vacxy,v%vacz,kvac,nv2,& + tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk) ! !---> generate a and b coeffficients ! - IF (noco%l_noco) THEN - DO jspin = 1,input%jspins + DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2) DO k = 1,lapw%nv(jspin) gz = sign*cell%bmat(3,3)*lapw%k3(k,jspin) i2 = map2(k,jspin) @@ -109,16 +103,6 @@ CONTAINS b(k,jspin) = c_1 * CMPLX(duz(i2,jspin), gz* uz(i2,jspin) ) ENDDO ENDDO - ELSE - DO k = 1,lapw%nv(jsp) - gz = sign*cell%bmat(3,3)*lapw%gvec(3,k,jsp) - i2 = map2(k,jsp) - th = gz*cell%z1 - c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk) - a(k,jsp) = - c_1 * CMPLX(dudz(i2,jsp), gz*udz(i2,jsp) ) - b(k,jsp) = c_1 * CMPLX(duz(i2,jsp), gz* uz(i2,jsp) ) - ENDDO - ENDIF !---> update hamiltonian and overlap matrices IF (jspin1==jspin2) THEN DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size @@ -153,30 +137,25 @@ CONTAINS !Diagonal term of Overlapp matrix, Hamiltonian later sij = CONJG(a(i,jspin2))*a(i,jspin2) + CONJG(b(i,jspin2))*b(i,jspin2)*ddnv(ik,jspin2) IF (hmat(1,1)%l_real) THEN - smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij) + smat(s2,s1)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij) ELSE - smat(s1,s2)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij + smat(s2,s1)%data_c(j,i0) = smat(s1,s2)%data_c(j,i0) + sij ENDIF ENDDO ENDIF !---> hamiltonian update - DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size + DO i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size i0=(i-1)/mpi%n_size+1 !local column index - ik = map2(i,jspin2) - DO j = 1,MIN(i,lapw%nv(jspin1)) - jk = map2(j,jspin1) - IF (jspin2>jspin1) THEN - hij = CONJG(CONJG(a(j,jspin1))* (tuuv(jk,ik)*a(i,jspin2) +tudv(jk,ik)*b(i,jspin2))& - + CONJG(b(j,jspin1))* (tddv(jk,ik)*b(i,jspin2) +tduv(jk,ik)*a(i,jspin2))) - ELSE - hij = CONJG(a(i,jspin2))* (tuuv(ik,jk)*a(j,jspin1) +tudv(ik,jk)*b(j,jspin1))& - + CONJG(b(i,jspin2))* (tddv(ik,jk)*b(j,jspin1) +tduv(ik,jk)*a(j,jspin1)) - ENDIF + ik = map2(i,jspin1) + DO j = 1,MERGE(i,lapw%nv(jspin2),jspin1==jspin2) + jk = map2(j,jspin2) + hij = CONJG(a(i,jspin1))* (tuuv(ik,jk)*a(j,jspin2) +tudv(ik,jk)*b(j,jspin2))& + + CONJG(b(i,jspin1))* (tddv(ik,jk)*b(j,jspin2) +tduv(ik,jk)*a(j,jspin2)) IF (hmat(1,1)%l_real) THEN - hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + REAL(hij) + hmat(s2,s1)%data_r(j,i0) = hmat(s2,s1)%data_r(j,i0) + REAL(hij) ELSE - hmat(s1,s2)%data_c(j,i0) = hmat(s1,s2)%data_c(j,i0) + hij + hmat(s2,s1)%data_c(j,i0) = hmat(s2,s1)%data_c(j,i0) + hij ENDIF ENDDO ENDDO diff --git a/eigen/vacfun.f90 b/eigen/vacfun.f90 index 4fad04570b39d16dc208f2012be2f744fc7dfe59..e83399645f3506acc95f56cda9e1015f4c4bdcb9 100644 --- a/eigen/vacfun.f90 +++ b/eigen/vacfun.f90 @@ -2,8 +2,8 @@ MODULE m_vacfun use m_juDFT CONTAINS SUBROUTINE vacfun(& - vacuum,DIMENSION,stars, jsp,input,noco,jsp1,jsp2,& - sym, cell,ivac,evac,bkpt, vxy,vz,kvac1,kvac2,nv2,& + vacuum,stars,input,noco,jspin1,jspin2,& + sym, cell,ivac,evac,bkpt, vxy,vz,kvac,nv2,& tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk) !********************************************************************* ! determines the necessary values and derivatives on the vacuum @@ -19,7 +19,6 @@ CONTAINS USE m_types IMPLICIT NONE - TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_input),INTENT(IN) :: input TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_noco),INTENT(IN) :: noco @@ -28,31 +27,35 @@ CONTAINS TYPE(t_cell),INTENT(IN) :: cell ! .. ! .. Scalar Arguments .. - INTEGER, INTENT (IN) :: jsp ,ivac,jsp1,jsp2 + INTEGER, INTENT (IN) :: ivac,jspin1,jspin2 REAL, INTENT (OUT) :: wronk ! .. ! .. Array Arguments .. - INTEGER, INTENT (IN) :: nv2(input%jspins) - INTEGER, INTENT (IN) :: kvac1(dimension%nv2d,input%jspins),kvac2(dimension%nv2d,input%jspins) - COMPLEX, INTENT (IN) :: vxy(vacuum%nmzxyd,stars%ng2-1) - COMPLEX, INTENT (OUT):: tddv(dimension%nv2d,dimension%nv2d),tduv(dimension%nv2d,dimension%nv2d) - COMPLEX, INTENT (OUT):: tudv(dimension%nv2d,dimension%nv2d),tuuv(dimension%nv2d,dimension%nv2d) - REAL, INTENT (IN) :: vz(vacuum%nmzd,2,4) ,evac(2,input%jspins) + INTEGER, INTENT (IN) :: nv2(:)!(input%jspins) + INTEGER, INTENT (IN) :: kvac(:,:,:)!(2,dimension%nv2d,input%jspins) + COMPLEX, INTENT (IN) :: vxy(:,:,:,:) !(vacuum%nmzxyd,stars%ng2-1,nvac,:) + COMPLEX, INTENT (OUT):: tddv(:,:),tduv(:,:)!(dimension%nv2d,dimension%nv2d) + COMPLEX, INTENT (OUT):: tudv(:,:),tuuv(:,:)!(dimension%nv2d,dimension%nv2d) + REAL,ALLOCATABLE,INTENT (IN) :: vz(:,:,:) !(vacuum%nmzd,2,4) , + REAL, INTENT (IN) :: evac(:,:)!(2,input%jspins) REAL, INTENT (IN) :: bkpt(3) - REAL, INTENT (OUT):: udz(dimension%nv2d,input%jspins),uz(dimension%nv2d,input%jspins) - REAL, INTENT (OUT):: dudz(dimension%nv2d,input%jspins),duz(dimension%nv2d,input%jspins) - REAL, INTENT (OUT):: ddnv(dimension%nv2d,input%jspins) + REAL, INTENT (OUT):: udz(:,:),uz(:,:)!(dimension%nv2d,input%jspins) + REAL, INTENT (OUT):: dudz(:,:),duz(:,:)!(dimension%nv2d,input%jspins) + REAL, INTENT (OUT):: ddnv(:,:)!(dimension%nv2d,input%jspins) ! .. ! .. Local Scalars .. - REAL ev,scale,xv,yv,vzero + REAL ev,scale,xv,yv,vzero,fac COMPLEX phase - INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin + INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin,ipot LOGICAL tail ! .. ! .. Local Arrays .. - REAL u(vacuum%nmzd,dimension%nv2d,input%jspins),ud(vacuum%nmzd,dimension%nv2d,input%jspins) + REAL u(vacuum%nmzd,size(duz,1),input%jspins),ud(vacuum%nmzd,size(duz,1),input%jspins) REAL v(3),x(vacuum%nmzd), qssbti(2,2) ! .. + fac=MERGE(1.0,-1.0,jspin1<=jspin2) + ipot=MERGE(jspin1,3,jspin1==jspin2) + tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0 udz=0.0;duz=0.0;ddnv=0.0;udz=0.;uz=0. tail = .true. @@ -60,15 +63,12 @@ CONTAINS !---> wronksian for the schrodinger equation given by an identity wronk = 2.0 !---> setup the spin-spiral q-vector - qssbti(1,1) = - noco%qss(1)/2 - qssbti(2,1) = - noco%qss(2)/2 - qssbti(1,2) = + noco%qss(1)/2 - qssbti(2,2) = + noco%qss(2)/2 + qssbti(1:2,1) = - noco%qss(1:2)/2 + qssbti(1:2,2) = + noco%qss(1:2)/2 !---> generate basis functions for each 2-d k+g - DO jspin = 1,input%jspins + DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2) DO ik = 1,nv2(jspin) - v(1) = bkpt(1) + kvac1(ik,jspin) + qssbti(1,jspin) - v(2) = bkpt(2) + kvac2(ik,jspin) + qssbti(2,jspin) + v(1:2) = bkpt(1:2) + kvac(:,ik,jspin) + qssbti(1:2,jspin) v(3) = 0.0 ev = evac(ivac,jspin) - 0.5*dot_product(v,matmul(v,cell%bbmat)) vzero = vz(vacuum%nmzd,ivac,jspin) @@ -80,19 +80,19 @@ CONTAINS !---> make sure the solutions satisfy the wronksian scale = wronk/ (udz(ik,jspin)*duz(ik,jspin)-& & dudz(ik,jspin)*uz(ik,jspin)) - udz(ik,jspin) = scale*udz(ik,jspin) + udz(ik,jspin) = scale*udz(ik,jspin) dudz(ik,jspin) = scale*dudz(ik,jspin) ddnv(ik,jspin) = scale*ddnv(ik,jspin) ud(:,ik,jspin) = scale*ud(:,ik,jspin) enddo ENDDO !---> set up the tuuv, etc. matrices - DO ik = 1,nv2(jsp1) - DO jk = 1,nv2(jsp2) + DO ik = 1,nv2(jspin1) + DO jk = 1,nv2(jspin2) !---> determine the warping component of the potential - i1 = kvac1(ik,jsp1) - kvac1(jk,jsp2) - i2 = kvac2(ik,jsp1) - kvac2(jk,jsp2) + i1 = kvac(1,ik,jspin1) - kvac(1,jk,jspin2) + i2 = kvac(2,ik,jspin1) - kvac(2,jk,jspin2) i3 = 0 ind3 = stars%ig(i1,i2,i3) IF (ind3.EQ.0) CYCLE @@ -114,98 +114,96 @@ CONTAINS !---> tuuv DO i = 1,vacuum%nmzxy - x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2)) + x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail) DO i = 1,vacuum%nmzxy - x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2)) + x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail) tuuv(ik,jk) = phase*cmplx(xv,yv) !---> tddv DO i = 1,vacuum%nmzxy - x(np1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2)) + x(np1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail) DO i = 1,vacuum%nmzxy - x(np1-i) =ud(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2)) + x(np1-i) =ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail) tddv(ik,jk) = phase*cmplx(xv,yv) !---> tudv DO i = 1,vacuum%nmzxy - x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2)) + x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail) DO i = 1,vacuum%nmzxy - x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2)) + x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail) tudv(ik,jk) = phase*cmplx(xv,yv) !---> tduv DO i = 1,vacuum%nmzxy - x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2)) + x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail) DO i = 1,vacuum%nmzxy - x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2)) + x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot)) enddo CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail) tduv(ik,jk) = phase*cmplx(xv,yv) ELSE - !---> diagonal (film muffin-tin) terms - IF (jsp1==jsp2) THEN - tuuv(ik,ik) = cmplx(evac(ivac,jsp1),0.0) - tddv(ik,ik) = cmplx(evac(ivac,jsp1)*ddnv(ik,jsp1),0.0) + IF (jspin1==jspin2) THEN + tuuv(ik,ik) = cmplx(evac(ivac,jspin1),0.0) + tddv(ik,ik) = cmplx(evac(ivac,jspin1)*ddnv(ik,jspin1),0.0) tudv(ik,ik) = cmplx(0.5,0.0) tduv(ik,ik) = cmplx(0.5,0.0) ELSE - !---> tuuv DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3) + x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail) DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4) + x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail) tuuv(ik,jk) = cmplx(xv,yv) !---> tddv DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3) + x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail) DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4) + x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail) tddv(ik,jk) = cmplx(xv,yv) !---> tudv DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3) + x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail) DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4) + x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail) tudv(ik,jk) = cmplx(xv,yv) !---> tduv DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3) + x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail) DO i = 1,vacuum%nmz - x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4) + x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4) ENDDO CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail) tduv(ik,jk) = cmplx(xv,yv) diff --git a/force/CMakeLists.txt b/force/CMakeLists.txt index 9af5f733e278542cd6b520404f13d42cd1d42e7e..3658f5c62b10603aa3aa44da1df26629baded3f9 100644 --- a/force/CMakeLists.txt +++ b/force/CMakeLists.txt @@ -1,9 +1,9 @@ set(fleur_F77 ${fleur_F77} -force/bfgs.f -force/bfgs0.f force/relax.F ) set(fleur_F90 ${fleur_F90} +force/bfgs0.f90 +force/bfgs.f90 force/force_b8.f90 force/force_0.f90 force/force_a12.f90 diff --git a/force/bfgs.f b/force/bfgs.f deleted file mode 100644 index 7baeaa6cbeef6dfdcba5b46c6bf3fd064c492d4e..0000000000000000000000000000000000000000 --- a/force/bfgs.f +++ /dev/null @@ -1,197 +0,0 @@ - MODULE m_bfgs - use m_juDFT - CONTAINS - SUBROUTINE bfgs( - > ntype,istep,istep0,force, - > zat,xa,thetad,epsdisp,epsforce,tote, - X xold,y,h,tau0, - < lconv) -!******************************************************************* -! relaxes the forces using the BFGS quasi-Newton method. -! input: -! istep = atomic step in this run -! istep0 = number of atomic steps in previous runs used -! in updating hessian -! -! output: -! lconv = logical true if forces are converged to tolerance -! given in epsforce -! -! the positions and forces from this step are added -! to file force.dat -!******************************************************************* - - IMPLICIT NONE -! .. -! .. Scalar Arguments .. - INTEGER, INTENT (IN) :: ntype,istep,istep0 - REAL, INTENT (IN) :: epsdisp,epsforce,tote - REAL, INTENT (INOUT) :: thetad,xa - LOGICAL, INTENT (OUT) :: lconv -! .. -! .. Array Arguments .. - REAL, INTENT (IN) :: force(3,ntype),zat(ntype) - REAL, INTENT (INOUT) :: tau0(3,ntype) - REAL, INTENT (INOUT) :: y(3*ntype),xold(3*ntype) - REAL, INTENT (INOUT) :: h(3*ntype,3*ntype) -! .. -! .. Local Scalars .. - INTEGER :: i,j,n,nn,ist,na - REAL :: py,yy,alpha,s,zatm,fmax,gamma,d2,dispmax -! .. -! .. Local Arrays .. - REAL :: f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype) -! - n=3*ntype - ist = istep+istep0 - -!--- > get positions and forces in correct form and output - DO na = 1,ntype - nn = 3*(na-1) - xnew(nn+1:nn+3) = tau0(:,na) - f(nn+1:nn+3) = force(:,na) - ENDDO - !Write new entry into forces.dat - OPEN(43,file ='forces.dat',status ='unknown',form='formatted' - $ ,position ='append') - WRITE (43,'(a,f20.10)') "energy =",tote - WRITE (43,'(3f16.9,3f14.9)') - > ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1) - CLOSE(43) -!--- > get maximum force - fmax = 0.0 - DO na = 1,ntype - nn = 3*(na-1) - fmax = MAX( fmax, (f(nn+1)**2+f(nn+2)**2+f(nn+3)**2) ) - ENDDO - fmax = SQRT(fmax) - WRITE (6,1000) istep,fmax - IF ( fmax if first step, go along gradient - - IF (ist==1) THEN -!--- > choose a reasonable first guess for scaling, but -!--- > limit displacement to a maximum of 0.25 a.u. -!--- > (may need to be changed for different systems) -!--- > this choice is based on a Debye temperature of 330K; -!--- > modify as needed (change thetad in param.8) - zatm = 0.0 - DO i = 1,ntype - zatm = MAX(zatm,zat(i)) - ENDDO - IF (ABS(xa)<1.0e-10) THEN - WRITE (6,*) 'WARNING, xa = 0.0 set to 2.0' - xa = 2.0 - ENDIF - IF (ABS(thetad)<1.0e-10) THEN - WRITE (6,*) 'WARNING, thetad = 0.0 set to 330.0' - thetad = 330.0 - ENDIF - alpha = (250.0/(zatm*xa))*((330./thetad)**2) - IF ( alpha*fmax*xa > 0.15 ) alpha = 0.25/(fmax*xa) - p(:) = alpha*f(:) -!--- > set h to identity matrix - h = 0.0 - DO j = 1,n - h(j,j) = 1.0 - ENDDO - - ELSE -!------------------------------------------------------------ -! ===> k-th step -! -!--- > determine p - p(:) = xnew(:)-xold(:) -!--- > update the change in gradients - y(:) = y(:)-f(:) -!--- > get necessary inner products and H|y> - py = dot_PRODUCT(p,y) - yy = 0.0 - DO i = 1,n - s = 0.0 - DO j = 1,n - s = s+y(j)*h(j,i) - ENDDO - v(i) = s - yy = yy+y(i)*s - ENDDO -!--- > check that update will leave h positive definite; -!--- > if not, then stop - IF (py<=0.0) THEN - WRITE (6,*) ' bfgs: < 0' - WRITE (6,*) ' check convergence of forces' - CALL juDFT_error("bfgs: <0",calledby="bfgs") - ELSE -!--- > update h - IF (ist==2) THEN - gamma = py/yy - ELSE - gamma = 1.0 - ENDIF - DO j = 1,n - DO i = 1,n - h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma - + + (1.+gamma*yy/py)*p(i)*p(j)/py - ENDDO - ENDDO -!--- > generate p - DO i = 1,n - s = 0.0 - DO j = 1,n - s = s+f(j)*h(j,i) - ENDDO - p(i) = s - ENDDO - ENDIF - ENDIF - -!------------------------------------------------------------- -!--- > put xold and y in the correct form for the next step - DO i = 1,n - xold(i) = xnew(i) - y(i) = f(i) - ENDDO - -!--- > if displacements are all less than epsdisp, then converged - dispmax = 0.0 - DO na = 1,ntype - nn = 3*(na-1) - d2 = p(nn+1)**2 + p(nn+2)**2 + p(nn+3)**2 - dispmax = MAX( dispmax, d2) - ENDDO - dispmax = xa*SQRT(dispmax) - IF (dispmax get new displacements - DO i = 1,n - xnew(i) = xold(i)+p(i) - ENDDO - DO na = 1,ntype - nn = 3*(na-1) - tau0(1,na) = xnew(nn+1) - tau0(2,na) = xnew(nn+2) - tau0(3,na) = xnew(nn+3) - ENDDO - - WRITE (6,'(1x/)') - WRITE (6,*) 'changes in p for step',ist - DO na = 1,ntype - nn = 3*(na-1) - WRITE (6,'(i5,6f12.6)') na,(p(nn+i),i = 1,3),(xnew(nn+j),j=1,3) - ENDDO - - END SUBROUTINE bfgs - END MODULE m_bfgs diff --git a/force/bfgs.f90 b/force/bfgs.f90 new file mode 100644 index 0000000000000000000000000000000000000000..95c994de7a5a18860295b94bd36480348a382967 --- /dev/null +++ b/force/bfgs.f90 @@ -0,0 +1,191 @@ +MODULE m_bfgs + +USE m_juDFT + +CONTAINS + +!******************************************************************* +! relaxes the forces using the BFGS quasi-Newton method. +! input: +! istep = atomic step in this run +! istep0 = number of atomic steps in previous runs used +! in updating hessian +! +! output: +! lconv = logical true if forces are converged to tolerance +! given in epsforce +! +! the positions and forces from this step are added +! to file force.dat +!******************************************************************* +SUBROUTINE bfgs(ntype,istep,istep0,force,zat,xa,thetad,epsdisp,epsforce,tote,& + xold,y,h,tau0,lconv) + +IMPLICIT NONE + + ! Scalar Arguments + INTEGER, INTENT (IN) :: ntype,istep,istep0 + REAL, INTENT (IN) :: epsdisp,epsforce,tote + REAL, INTENT (INOUT) :: thetad,xa + LOGICAL, INTENT (OUT) :: lconv + + ! Array Arguments + REAL, INTENT (IN) :: force(3,ntype),zat(ntype) + REAL, INTENT (INOUT) :: tau0(3,ntype) + REAL, INTENT (INOUT) :: y(3*ntype),xold(3*ntype) + REAL, INTENT (INOUT) :: h(3*ntype,3*ntype) + + ! Local Scalars + INTEGER :: i,j,n,nn,ist,na + REAL :: py,yy,alpha,s,zatm,fmax,gamma,d2,dispmax + + ! Local Arrays + REAL :: f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype) + + n=3*ntype + ist = istep+istep0 + + ! get positions and forces in correct form and output + DO na = 1,ntype + nn = 3*(na-1) + xnew(nn+1:nn+3) = tau0(:,na) + f(nn+1:nn+3) = force(:,na) + END DO + + ! Write new entry into forces.dat + OPEN(43,file ='forces.dat',status ='unknown',form='formatted',position ='append') + WRITE (43,'(a,f20.10)') "energy =",tote + WRITE (43,'(3f16.9,3f14.9)') ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1) + CLOSE(43) + + ! get maximum force + fmax = 0.0 + DO na = 1,ntype + nn = 3*(na-1) + fmax = MAX( fmax, (f(nn+1)**2+f(nn+2)**2+f(nn+3)**2) ) + END DO + fmax = SQRT(fmax) + WRITE (6,1000) istep,fmax + IF (fmax < epsforce) THEN + lconv = .TRUE. + RETURN + ELSE + lconv = .FALSE. + END IF + 1000 FORMAT (1x/,' atomic step',i4,': maximum force =',1p,e14.6,' hartrees/a.u.') + + ! if first step, go along gradient + IF (ist==1) THEN + ! choose a reasonable first guess for scaling, but + ! limit displacement to a maximum of 0.25 a.u. + ! (may need to be changed for different systems) + ! this choice is based on a Debye temperature of 330K; + ! modify as needed (change thetad in param.8) + zatm = 0.0 + DO i = 1,ntype + zatm = MAX(zatm,zat(i)) + END DO + IF (ABS(xa)<1.0e-10) THEN + WRITE (6,*) 'WARNING, xa = 0.0 set to 2.0' + xa = 2.0 + END IF + IF (ABS(thetad)<1.0e-10) THEN + WRITE (6,*) 'WARNING, thetad = 0.0 set to 330.0' + thetad = 330.0 + END IF + alpha = (250.0/(zatm*xa))*((330./thetad)**2) + IF ( alpha*fmax*xa > 0.15 ) alpha = 0.25/(fmax*xa) + p(:) = alpha*f(:) + ! set h to identity matrix + h = 0.0 + DO j = 1,n + h(j,j) = 1.0 + END DO + ELSE + ! k-th step + + ! determine p + p(:) = xnew(:) - xold(:) + ! update the change in gradients + y(:) = y(:) - f(:) + ! get necessary inner products and H|y> + py = dot_PRODUCT(p,y) + yy = 0.0 + DO i = 1,n + s = 0.0 + DO j = 1,n + s = s+y(j)*h(j,i) + END DO + v(i) = s + yy = yy+y(i)*s + END DO + ! check that update will leave h positive definite; + ! if not, then stop + IF (py <= 0.0) THEN + WRITE (6,*) ' bfgs: < 0' + WRITE (6,*) ' check convergence of forces' + CALL juDFT_error("bfgs: <0",calledby="bfgs") + ELSE + ! update h + IF (ist==2) THEN + gamma = py/yy + ELSE + gamma = 1.0 + ENDIF + DO j = 1,n + DO i = 1,n + h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma + (1.+gamma*yy/py)*p(i)*p(j)/py + END DO + END DO + ! generate p + DO i = 1,n + s = 0.0 + DO j = 1,n + s = s+f(j)*h(j,i) + END DO + p(i) = s + END DO + END IF + END IF + + ! put xold and y in the correct form for the next step + DO i = 1, n + xold(i) = xnew(i) + y(i) = f(i) + END DO + + ! if displacements are all less than epsdisp, then converged + dispmax = 0.0 + DO na = 1,ntype + nn = 3*(na-1) + d2 = p(nn+1)**2 + p(nn+2)**2 + p(nn+3)**2 + dispmax = MAX( dispmax, d2) + END DO + dispmax = xa*SQRT(dispmax) + IF (dispmax ntype, - < istep0,xold,y,h) -c -c******************************************************************* -c checks whether a file exists with positions and forces to -c use to relax structure using the BFGS quasi-Newton method. -c istep0 is the number of steps included on file. -c h,y, and xold will be saved correctly for further use. -c******************************************************************* -c - IMPLICIT NONE -C .. -C .. Scalar Arguments .. - INTEGER, INTENT (IN) :: ntype - INTEGER, INTENT (OUT):: istep0 -C .. -C .. Array Arguments .. - REAL, INTENT (OUT):: y(3*ntype),xold(3*ntype) - REAL, INTENT (OUT):: h(3*ntype,3*ntype) -C .. -C .. Local Scalars .. - INTEGER n,is,i,j - REAL yy,py,s,gamma - CHARACTER :: testchar - INTEGER :: maxhist - LOGICAL :: l_skip -! .. -C .. Local Arrays .. - REAL f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype) -c -c---> open the file with forces -c - OPEN (43,file='forces.dat',status='unknown',form='formatted') - REWIND 43 - - n=3*ntype - istep0=0 - - !The might be a line like - !maxhist=3 - !at the beginning of the file - READ(43,'(a)',end=100) testchar - REWIND 43 - IF (testchar =='m') THEN - READ(43,"(8x,i5)") maxhist - DO - READ(43,*,END = 99) !read the line containing the energy - istep0 = istep0+1 - DO j = 1,ntype - READ(43,*,err= 100, end=100) - ENDDO - ENDDO - 99 REWIND 43 - read(43,*) !skip maxhist line - !skip all but last maxhist entries - DO i = 1, MAX(istep0-maxhist,0) - read(43,*,end=100) - DO j = 1,ntype - READ(43,*,err= 100, end=100) - ENDDO - ENDDO - ENDIF - istep0 = 0 - - DO ! read from file until end of data - l_skip = .FALSE. - READ(43,'(33x,l1)',END = 101,err = 101) l_skip - 101 READ (43,'(3f16.9,3f14.9)',END = 100,err = 100) - > ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1) - IF (l_skip) CYCLE !skip this entry - - istep0 = istep0+1 - - IF (istep0.EQ.1) THEN -!------------------------------------------------------------ -! ===> first step -! -!--- > set h to identity matrix - h = 0.0 - DO j = 1,n - h(j,j) = 1.0 - ENDDO - -c------------------------------------------------------------ - ELSE -! ===> k-th step - -!--- > determine p and shift x as needed - p(:) = xnew(:)-xold(:) -!--- > update the change in gradients - y(:) = y(:)-f(:) -!--- > get necessary inner products and H|y> - py = dot_product(p,y) - yy = 0.0 - DO i = 1,n - s = 0.0 - DO j = 1,n - s = s+y(j)*h(j,i) - ENDDO - v(i) = s - yy = yy+y(i)*s - ENDDO -!--- > check that update will leave h positive definite, -!--- > if not, restart with diagonal matrix - IF (py<=0.0) THEN - WRITE (6,*) 'bfgs0: < 0, restart istep0 =',istep0 - h = 0.0 - DO j = 1,n - h(j,j) = 1.0 - ENDDO - istep0 = 1 - ELSE -!--- > for second step, use oren-spedicato scaling -!--- > for standard scaling, use gamma = 1. - IF (istep0 == 2) THEN - gamma = py/yy - ELSE - gamma = 1.0 - ENDIF -!--- > update h - DO j = 1,n - DO i = 1,n - h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma - + + (1.+gamma*yy/py)*p(i)*p(j)/py - ENDDO - ENDDO - ENDIF - ENDIF - -!------------------------------------------------------------ -!--- > put y in the correct form for the next step - DO i = 1,n - xold(i) = xnew(i) - y(i) = f(i) - ENDDO - - ENDDO - - 100 CLOSE(43) - IF (istep0==0) THEN - WRITE (6,1000) - ELSE - WRITE (6,1001) istep0 - ENDIF - 1000 FORMAT(' bfgs0: No previous atomic relaxation data found'/) - 1001 FORMAT(' bfgs0: Approximate hessian constructed using',i4, - + ' steps'/) - - END SUBROUTINE bfgs0 - END MODULE m_bfgs0 diff --git a/force/bfgs0.f90 b/force/bfgs0.f90 new file mode 100644 index 0000000000000000000000000000000000000000..dbdd5d6b15a45a8fe57f782c1157f7c4bece0e1c --- /dev/null +++ b/force/bfgs0.f90 @@ -0,0 +1,148 @@ +MODULE m_bfgs0 + +CONTAINS + +!******************************************************************* +! checks whether a file exists with positions and forces to +! use to relax structure using the BFGS quasi-Newton method. +! istep0 is the number of steps included on file. +! h,y, and xold will be saved correctly for further use. +!******************************************************************* +SUBROUTINE bfgs0(ntype,istep0,xold,y,h) + + IMPLICIT NONE + + ! Scalar Arguments + INTEGER, INTENT (IN) :: ntype + INTEGER, INTENT (OUT):: istep0 + + ! Array Arguments + REAL, INTENT (OUT):: y(3*ntype),xold(3*ntype) + REAL, INTENT (OUT):: h(3*ntype,3*ntype) + + ! Local Scalars + INTEGER n,is,i,j + REAL yy,py,s,gamma + CHARACTER :: testchar + INTEGER :: maxhist + LOGICAL :: l_skip + + ! Local Arrays + REAL f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype) + + ! open the file with forces + + OPEN (43,file='forces.dat',status='unknown',form='formatted') + REWIND 43 + + n=3*ntype + istep0=0 + + ! There might be a line like + ! maxhist=3 + ! at the beginning of the file + READ(43,'(a)',end=100) testchar + REWIND 43 + IF (testchar =='m') THEN + READ(43,"(8x,i5)") maxhist + DO + READ(43,*,END = 99) !read the line containing the energy + istep0 = istep0+1 + DO j = 1,ntype + READ(43,*,err= 100, end=100) + END DO + END DO + 99 REWIND 43 + READ(43,*) !skip maxhist line + !skip all but last maxhist entries + DO i = 1, MAX(istep0-maxhist,0) + READ(43,*,end=100) + DO j = 1,ntype + READ(43,*,err= 100, end=100) + END DO + END DO + END IF + istep0 = 0 + + DO ! read from file until end of data + l_skip = .FALSE. + READ(43,'(33x,l1)',END = 101,err = 101) l_skip + 101 READ (43,'(3f16.9,3f14.9)',END = 100,err = 100) ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1) + IF (l_skip) CYCLE !skip this entry + + istep0 = istep0+1 + + IF (istep0.EQ.1) THEN + ! first step + + ! set h to identity matrix + h = 0.0 + DO j = 1,n + h(j,j) = 1.0 + END DO + + ELSE + ! k-th step + + ! determine p and shift x as needed + p(:) = xnew(:) - xold(:) + ! update the change in gradients + y(:) = y(:) - f(:) + ! get necessary inner products and H|y> + py = dot_product(p,y) + yy = 0.0 + DO i = 1,n + s = 0.0 + DO j = 1,n + s = s+y(j)*h(j,i) + END DO + v(i) = s + yy = yy+y(i)*s + END DO + ! check that update will leave h positive definite, if not, restart with diagonal matrix + IF (py<=0.0) THEN + WRITE (6,*) 'bfgs0: < 0, restart istep0 =',istep0 + h = 0.0 + DO j = 1,n + h(j,j) = 1.0 + END DO + istep0 = 1 + ELSE + ! for second step, use oren-spedicato scaling + ! for standard scaling, use gamma = 1. + IF (istep0 == 2) THEN + gamma = py/yy + ELSE + gamma = 1.0 + END IF + ! update h + DO j = 1,n + DO i = 1,n + h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma + (1.+gamma*yy/py)*p(i)*p(j)/py + END DO + END DO + END IF + END IF + + ! put y in the correct form for the next step + DO i = 1,n + xold(i) = xnew(i) + y(i) = f(i) + END DO + + END DO + + 100 CLOSE(43) + + IF (istep0==0) THEN + WRITE (6,1000) + ELSE + WRITE (6,1001) istep0 + END IF + + 1000 FORMAT(' bfgs0: No previous atomic relaxation data found'/) + 1001 FORMAT(' bfgs0: Approximate hessian constructed using',i4,' steps'/) + + END SUBROUTINE bfgs0 + +END MODULE m_bfgs0 diff --git a/global/chkmt.f90 b/global/chkmt.f90 index 792d87581c43a3bdbc502c86b17333c531c67d03..ea1247749c4caa1c7dd620bcc3f304d5e44667aa 100644 --- a/global/chkmt.f90 +++ b/global/chkmt.f90 @@ -391,7 +391,11 @@ jri11 = NINT(220*rmt1(n)) END IF jri11 = NINT(jri11*0.5) * 2 + 1 - dx11 = LOG(3200*atoms%nz(n)*rmt1(n))/(jri11-1) + IF (atoms%nz(n) > 0) THEN + dx11 = LOG(3200*atoms%nz(n)*rmt1(n))/(jri11-1) + ELSE + dx11 = LOG(3200*rmt1(n))/(jri11-1) + ENDIF rkm = MAX(rkm, lmax11/rmt1(n)) WRITE (6,'(a3,i3,2i5,2f10.6)') noel(n),atoms%nz(n),lmax11,jri11,rmt1(n),dx11 dx1(n) = dx11 diff --git a/hybrid/coulombmatrix.F90 b/hybrid/coulombmatrix.F90 index 02af57c4147e0c75b3bd1800ac9d2f10cda51937..c74a2016fac051f71b51303944a743ef786689ed 100644 --- a/hybrid/coulombmatrix.F90 +++ b/hybrid/coulombmatrix.F90 @@ -2213,14 +2213,14 @@ CONTAINS ! Ensure numerical stability. If both formulas are not sufficiently stable, the program stops. IF(r1.GT.r2) THEN IF(r1.LT.1d-6 .AND. l_warn) THEN - WRITE(6,'(A,E10.5,A,E10.5,A)') 'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:',r1,'(',r2,')' + WRITE(6,'(A,E12.5,A,E12.5,A)') 'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:',r1,'(',r2,')' WRITE(6,'(A,2F15.10,I4)') ' Current qnorms and atom type:', q1,q2,itype l_warned = .TRUE. END IF sphbessel_integral = s**3 / dq * da ELSE IF(r2.LT.1d-6 .AND. l_warn) THEN - WRITE(6,'(A,E10.5,A,E10.5,A)') 'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:',r2,'(',r1,')' + WRITE(6,'(A,E13.5,A,E13.5,A)') 'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:',r2,'(',r1,')' WRITE(6,'(A,2F15.10,I4)') ' Current qnorms and atom type:', q1,q2,itype l_warned = .TRUE. END IF diff --git a/init/mapatom.F90 b/init/mapatom.F90 index 15ae26f4d8ae60cd5a9f73dd3aa677cac9538cf3..ba61d6b6c1544cde4f7799a0a71c9b1de072fc74 100644 --- a/init/mapatom.F90 +++ b/init/mapatom.F90 @@ -136,7 +136,7 @@ CALL juDFT_error("mapatom",calledby="mapatom") ENDIF WRITE (6,FMT=8010) nat1,na,atoms%ngopr(na) - 8010 FORMAT (5x,'atom',i3,' can be mapped into atom',i3,& + 8010 FORMAT (5x,'atom',i5,' can be mapped into atom',i5,& & ' through group operation',i4) ! ! end of equivalent atoms diff --git a/init/postprocessInput.F90 b/init/postprocessInput.F90 index 139588c5931f07a3defe433124489e41c05826da..ac9a4b90e5a400237c9282d96602902233ce3595 100644 --- a/init/postprocessInput.F90 +++ b/init/postprocessInput.F90 @@ -530,7 +530,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts ! The following call to inpeig should not be required. ! CALL inpeig(atoms,cell,input,oneD%odd%d1,kpts,enpara) END IF - + kpts%wtkpt=kpts%wtkpt/sum(kpts%wtkpt) !Normalize k-point weight CALL prp_qfft(stars,cell,noco,input) diff --git a/io/writeBasis.F90 b/io/writeBasis.F90 index 464aadee9bb383eb08786ae39f5b05d7d9bc7e3d..8d54798a712b6e1909ef567e89076d89f4afd6e9 100644 --- a/io/writeBasis.F90 +++ b/io/writeBasis.F90 @@ -457,10 +457,10 @@ SUBROUTINE writeBasis(input,noco,kpts,atoms,sym,cell,enpara,vTot,vCoul,vx,mpi,DI IF(abs(bk(2)).LT.1e-7) bk(2) = abs(bk(2)) IF(abs(bk(3)).LT.1e-7) bk(3) = abs(bk(3)) !write(kpt_name , '(2a,i0)') TRIM(ADJUSTL(jsp_name)),'/kpt_',nk - write(kpt_name , '(a,f12.10,a,f12.10,a,f12.10)') '/kpt_',bk(1),',',bk(2),',',bk(3) - CALL h5lexists_f(jspGroupID, TRIM(ADJUSTL(kpt_name)), link_exists, hdfError) + write(kpt_name , '(2a,f12.10,a,f12.10,a,f12.10)') TRIM(ADJUSTL(jsp_name)),'/kpt_',bk(1),',',bk(2),',',bk(3) + CALL h5lexists_f(fileID, TRIM(ADJUSTL(kpt_name)), link_exists, hdfError) IF (link_exists) CYCLE - CALL h5gcreate_f(jspGroupID, TRIM(ADJUSTL(kpt_name)), kptGroupID, hdfError) + CALL h5gcreate_f(fileID, TRIM(ADJUSTL(kpt_name)), kptGroupID, hdfError) !--------------------abcoff, zmat, eig output here------------------- !,results%neig(nk,jsp),results%eig(:,nk,jsp) numbands=results%neig(nk,jsp) diff --git a/io/xmlOutput.F90 b/io/xmlOutput.F90 index 6a6ddb444a05064ed4e3ed619b4bf4d6450c1f05..c3491a1bb2481210b02e2841fbc0086f0f89fce0 100644 --- a/io/xmlOutput.F90 +++ b/io/xmlOutput.F90 @@ -22,6 +22,7 @@ MODULE m_xmlOutput INTEGER, SAVE :: currentElementIndex INTEGER, SAVE :: maxNumElements INTEGER, SAVE :: xmlOutputUnit + LOGICAL, SAVE :: xmlOpened = .FALSE. CHARACTER(LEN= 40), ALLOCATABLE :: elementList(:) PUBLIC startXMLOutput, endXMLOutput @@ -86,6 +87,7 @@ MODULE m_xmlOutput ELSE OPEN (xmlOutputUnit,file='out.xml',form='formatted',status='unknown') ENDIF + xmlOpened = .TRUE. WRITE (xmlOutputUnit,'(a)') '' WRITE (xmlOutputUnit,'(a)') '' CALL openXMLElement('programVersion',(/'version'/),(/version_const/)) @@ -150,6 +152,7 @@ MODULE m_xmlOutput CHARACTER(LEN=10) :: dateString CHARACTER(LEN=10) :: timeString + IF (.NOT.xmlOpened) RETURN DO WHILE (currentElementIndex.NE.0) CALL closeXMLElement(elementList(currentElementIndex)) END DO diff --git a/juDFT/time.F90 b/juDFT/time.F90 index e7ba4cf5e49dbe14fcead3c6e54692eb3d7efec8..c7f3ff9578af8c230c4dfa34cc1f6ed05619fa9c 100644 --- a/juDFT/time.F90 +++ b/juDFT/time.F90 @@ -245,7 +245,7 @@ CONTAINS ENDIF IF (time1) THEN - WRITE(fid,"(a,T7,a,T46,a,T85,a,i0,a,f9.0,a,f9.0,a)") TRIM(indentstring),TRIM(timername),TRIM(timestring(time,globaltimer%time,level)),"(",timer%no_calls,"calls:",timer%mintime,"sec -",timer%maxtime,"sec)" + WRITE(fid,"(a,T7,a,T46,a,T85,a,i0,a,f9.3,a,f9.3,a)") TRIM(indentstring),TRIM(timername),TRIM(timestring(time,globaltimer%time,level)),"% (",timer%no_calls,"calls:",timer%mintime,"sec -",timer%maxtime,"sec)" ELSE WRITE(fid,"(a,T7,a,T46,a)") TRIM(indentstring),TRIM(timername),timestring(time,globaltimer%time,level) END IF @@ -304,6 +304,7 @@ CONTAINS OPEN(2,FILE ='juDFT_times',STATUS="replace") ENDIF globaltimer%time=cputime()-globaltimer%starttime + globaltimer%starttime=cputime() WRITE(fn,"('Total execution time: ',i0,'sec')") INT(globaltimer%time) CALL add_usage_data("Runtime",globaltimer%time) CALL priv_writetimes_longest(globaltimer,fid=fn) diff --git a/optional/stden.f90 b/optional/stden.f90 index 921d313535fe64530ad8be1ae6c6cb5aa4221303..e849ef118bcd02f65bae19013a01e3a375c0a57a 100644 --- a/optional/stden.f90 +++ b/optional/stden.f90 @@ -48,7 +48,7 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,& TYPE(t_xcpot_inbuild) :: xcpot_dummy ! Local Scalars - REAL d,del,fix,h,r,rnot,z,bm,qdel + REAL d,del,fix,h,r,rnot,z,bm,qdel,va REAL denz1(1,1),vacxpot(1,1),vacpot(1,1) INTEGER i,ivac,iza,j,jr,k,n,n1,ispin INTEGER nw,ilo,natot,nat @@ -142,9 +142,10 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,& !---> new atom rnot = atoms%rmsh(1,n) IF (z.LT.1.0) THEN + va = max(z,1.e-8)/(input%jspins*sfp_const*atoms%volmts(n)) DO ispin = 1, input%jspins DO i = 1,jrc(n) ! dimension%msh - rh(i,n,ispin) = 1.e-10 + rh(i,n,ispin) = va/rat(i,n)**2 END DO END DO ELSE diff --git a/types/types_potden.F90 b/types/types_potden.F90 index 2884345fee1c7f0ca85c900166d5c2bc0d7b8f5c..ba6042b53cf4e343d113ea90485ae033b0362df8 100644 --- a/types/types_potden.F90 +++ b/types/types_potden.F90 @@ -10,6 +10,7 @@ MODULE m_types_potden INTEGER :: iter INTEGER :: potdenType COMPLEX,ALLOCATABLE :: pw(:,:),pw_w(:,:) + ! mt(radial_grid, sphhar, atom, spin) REAL,ALLOCATABLE :: mt(:,:,:,:) REAL,ALLOCATABLE :: vacz(:,:,:) COMPLEX,ALLOCATABLE :: vacxy(:,:,:,:) diff --git a/vgen/mt_tofrom_grid.F90 b/vgen/mt_tofrom_grid.F90 index a9054442824866adef74977559c19a59126812b0..9c1046d65fee4771871479e99aec79d29e810929 100644 --- a/vgen/mt_tofrom_grid.F90 +++ b/vgen/mt_tofrom_grid.F90 @@ -4,186 +4,185 @@ ! of the MIT license as expressed in the LICENSE file in more detail. !-------------------------------------------------------------------------------- MODULE m_mt_tofrom_grid - USE m_types - PRIVATE - REAL,PARAMETER :: d_15 = 1.e-15 - INTEGER,PARAMETER :: ndvgrd=6 ! this should be consistent across GGA derivative routines - REAL, ALLOCATABLE :: ylh(:,:,:),ylht(:,:,:),ylhtt(:,:,:) - REAL, ALLOCATABLE :: ylhf(:,:,:),ylhff(:,:,:),ylhtf(:,:,:) - REAL, ALLOCATABLE :: wt(:),rx(:,:),thet(:) - PUBLIC :: init_mt_grid,mt_to_grid,mt_from_grid,finish_mt_grid + USE m_types + PRIVATE + REAL, PARAMETER :: d_15 = 1.e-15 + INTEGER, PARAMETER :: ndvgrd = 6 ! this should be consistent across GGA derivative routines + REAL, ALLOCATABLE :: ylh(:, :, :), ylht(:, :, :), ylhtt(:, :, :) + REAL, ALLOCATABLE :: ylhf(:, :, :), ylhff(:, :, :), ylhtf(:, :, :) + REAL, ALLOCATABLE :: wt(:), rx(:, :), thet(:) + PUBLIC :: init_mt_grid, mt_to_grid, mt_from_grid, finish_mt_grid CONTAINS - SUBROUTINE init_mt_grid(nsp,jspins,atoms,sphhar,xcpot,sym) - USE m_gaussp - USE m_lhglptg - USE m_lhglpts - IMPLICIT NONE - INTEGER,INTENT(IN) :: nsp,jspins - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_sphhar),INTENT(IN) :: sphhar - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_sym),INTENT(IN) :: sym - - ! generate nspd points on a sherical shell with radius 1.0 - ! angular mesh equidistant in phi, - ! theta are zeros of the legendre polynomials - ALLOCATE(wt(nsp),rx(3,nsp),thet(nsp)) - CALL gaussp(atoms%lmaxd, rx,wt) - ! generate the lattice harmonics on the angular mesh - ALLOCATE ( ylh(nsp,0:sphhar%nlhd,sphhar%ntypsd)) - IF (xcpot%needs_grad()) THEN - ALLOCATE(ylht,MOLD=ylh ) - ALLOCATE(ylhtt,MOLD=ylh ) - ALLOCATE(ylhf,MOLD=ylh ) - ALLOCATE(ylhff,MOLD=ylh ) - ALLOCATE(ylhtf,MOLD=ylh ) - - CALL lhglptg(sphhar,atoms,rx,nsp,xcpot,sym,& - ylh,thet,ylht,ylhtt,ylhf,ylhff,ylhtf) - ELSE - CALL lhglpts( sphhar,atoms, rx,nsp, sym, ylh) - END IF - END SUBROUTINE init_mt_grid - - SUBROUTINE mt_to_grid(xcpot,jspins,atoms,sphhar,den_mt,nsp,n,grad,ch) - USE m_grdchlh - USE m_mkgylm - IMPLICIT NONE - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_sphhar),INTENT(IN) :: sphhar - REAL,INTENT(IN) :: den_mt(:,0:,:) - INTEGER,INTENT(IN) :: n,jspins,nsp - REAL,INTENT(OUT),OPTIONAL :: ch(:,:) - TYPE(t_gradients),INTENT(INOUT):: grad - - REAL, ALLOCATABLE :: chlh(:,:,:),chlhdr(:,:,:),chlhdrr(:,:,:) - REAL, ALLOCATABLE :: chdr(:,:),chdt(:,:),chdf(:,:),ch_tmp(:,:) - REAL, ALLOCATABLE :: chdrr(:,:),chdtt(:,:),chdff(:,:),chdtf(:,:) - REAL, ALLOCATABLE :: chdrt(:,:),chdrf(:,:) - INTEGER:: nd,lh,js,jr,kt,k - - nd = atoms%ntypsy(SUM(atoms%neq(:n-1))+1) - - ALLOCATE ( chlh(atoms%jmtd,0:sphhar%nlhd,jspins)) - ALLOCATE ( ch_tmp(nsp,jspins) ) - IF (xcpot%needs_grad()) THEN - ALLOCATE(chdr(nsp,jspins),chdt(nsp,jspins),chdf(nsp,jspins),chdrr(nsp,jspins),& - chdtt(nsp,jspins),chdff(nsp,jspins),chdtf(nsp,jspins),chdrt(nsp,jspins),& - chdrf(nsp,jspins) ) - ALLOCATE (chlhdr(atoms%jmtd,0:sphhar%nlhd,jspins)) - ALLOCATE (chlhdrr(atoms%jmtd,0:sphhar%nlhd,jspins)) - ENDIF - - DO lh = 0,sphhar%nlh(nd) - - ! calculates gradients of radial charge densities of l=> 0. - ! rho*ylh/r**2 is charge density. chlh=rho/r**2. - ! charge density=sum(chlh*ylh). - ! chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr. - - DO js = 1,jspins - DO jr = 1,atoms%jri(n) - chlh(jr,lh,js) = den_mt(jr,lh,js)/(atoms%rmsh(jr,n)*atoms%rmsh(jr,n)) - ENDDO - IF (xcpot%needs_grad()) CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),atoms%rmsh(1,n),& - chlh(1,lh,js),ndvgrd, chlhdr(1,lh,js),chlhdrr(1,lh,js)) - - ENDDO ! js - ENDDO ! lh - - kt=0 - DO jr = 1,atoms%jri(n) - ch_tmp(:,:) = 0.0 ! charge density (on extended grid for all jr) - ! following are at points on jr-th sphere. - ! generate the densities on an angular mesh - DO js = 1,jspins - DO lh = 0,sphhar%nlh(nd) - DO k = 1,nsp - ch_tmp(k,js) = ch_tmp(k,js) + ylh(k,lh,nd)*chlh(jr,lh,js) - ENDDO - ENDDO - ENDDO - IF (xcpot%needs_grad()) THEN - chdr(:,:) = 0.0 ! d(ch)/dr - chdt(:,:) = 0.0 ! d(ch)/dtheta - chdf(:,:) = 0.0 ! d(ch)/dfai - chdrr(:,:) = 0.0 ! dd(ch)/drr - chdtt(:,:) = 0.0 ! dd(ch)/dtt - chdff(:,:) = 0.0 ! dd(ch)/dff - chdtf(:,:) = 0.0 ! dd(ch)/dtf - chdrt(:,:) = 0.0 ! d(d(ch)/dr)dt - chdrf(:,:) = 0.0 ! d(d(ch)/dr)df - ! generate the derivatives on an angular mesh - DO js = 1,jspins - DO lh = 0,sphhar%nlh(nd) - ! - DO k = 1,nsp - chdr(k,js) =chdr(k,js)+ ylh(k,lh,nd)*chlhdr(jr,lh,js) - chdrr(k,js)=chdrr(k,js)+ylh(k,lh,nd)*chlhdrr(jr,lh,js) - ENDDO - - DO k = 1,nsp - chdrt(k,js)=chdrt(k,js)+ylht(k,lh,nd)*chlhdr(jr,lh,js) - chdrf(k,js)=chdrf(k,js)+ylhf(k,lh,nd)*chlhdr(jr,lh,js) - chdt(k,js) =chdt(k,js) +ylht(k,lh,nd)*chlh(jr,lh,js) - chdf(k,js) =chdf(k,js) +ylhf(k,lh,nd)*chlh(jr,lh,js) - chdtt(k,js)=chdtt(k,js)+ylhtt(k,lh,nd)*chlh(jr,lh,js) - chdff(k,js)=chdff(k,js)+ylhff(k,lh,nd)*chlh(jr,lh,js) - chdtf(k,js)=chdtf(k,js)+ylhtf(k,lh,nd)*chlh(jr,lh,js) - ENDDO - ENDDO ! lh - ENDDO ! js - - - CALL mkgylm(jspins,atoms%rmsh(jr,n),thet,nsp,& - ch_tmp,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,grad,kt) - ENDIF - !Set charge to minimum value - IF (PRESENT(ch)) THEN - WHERE(ABS(ch_tmp) < d_15) ch_tmp = d_15 - ch(kt+1:kt+nsp,:) = ch_tmp(:nsp,:) - ENDIF - kt=kt+nsp - END DO - - END SUBROUTINE mt_to_grid - - - SUBROUTINE mt_from_grid(atoms,sphhar,nsp,n,jspins,v_in,vr) - IMPLICIT NONE - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_sphhar),INTENT(IN):: sphhar - INTEGER,INTENT(IN) :: nsp,jspins,n - REAL,INTENT(IN) :: v_in(:,:) - REAL,INTENT(INOUT) :: vr(:,0:,:) - - REAL :: vpot(nsp),vlh - INTEGER :: js,kt,lh,jr,nd - nd=atoms%ntypsy(SUM(atoms%neq(:n-1))+1) - DO js = 1,jspins - ! - kt=0 - DO jr=1,atoms%jri(n) - vpot=v_in(kt+1:kt+nsp,js)*wt(:)! multiplicate v_in with the weights of the k-points - - DO lh = 0,sphhar%nlh(nd) - ! - ! ---> determine the corresponding potential number - !c through gauss integration - ! - vlh=dot_PRODUCT(vpot(:),ylh(:nsp,lh,nd)) - vr(jr,lh,js) = vr(jr,lh,js) + vlh - ENDDO ! lh - kt=kt+nsp - ENDDO ! jr - ENDDO - - END SUBROUTINE mt_from_grid - - SUBROUTINE finish_mt_grid() - DEALLOCATE(ylh,wt,rx,thet) - IF (ALLOCATED(ylht)) DEALLOCATE(ylht,ylhtt,ylhf,ylhff,ylhtf) - END SUBROUTINE finish_mt_grid + SUBROUTINE init_mt_grid(nsp, jspins, atoms, sphhar, xcpot, sym) + USE m_gaussp + USE m_lhglptg + USE m_lhglpts + IMPLICIT NONE + INTEGER, INTENT(IN) :: nsp, jspins + TYPE(t_atoms), INTENT(IN) :: atoms + TYPE(t_sphhar), INTENT(IN) :: sphhar + CLASS(t_xcpot), INTENT(IN) :: xcpot + TYPE(t_sym), INTENT(IN) :: sym + + ! generate nspd points on a sherical shell with radius 1.0 + ! angular mesh equidistant in phi, + ! theta are zeros of the legendre polynomials + ALLOCATE (wt(nsp), rx(3, nsp), thet(nsp)) + CALL gaussp(atoms%lmaxd, rx, wt) + ! generate the lattice harmonics on the angular mesh + ALLOCATE (ylh(nsp, 0:sphhar%nlhd, sphhar%ntypsd)) + IF (xcpot%needs_grad()) THEN + ALLOCATE (ylht, MOLD=ylh) + ALLOCATE (ylhtt, MOLD=ylh) + ALLOCATE (ylhf, MOLD=ylh) + ALLOCATE (ylhff, MOLD=ylh) + ALLOCATE (ylhtf, MOLD=ylh) + + CALL lhglptg(sphhar, atoms, rx, nsp, xcpot, sym, & + ylh, thet, ylht, ylhtt, ylhf, ylhff, ylhtf) + ELSE + CALL lhglpts(sphhar, atoms, rx, nsp, sym, ylh) + END IF + END SUBROUTINE init_mt_grid + + SUBROUTINE mt_to_grid(xcpot, jspins, atoms, sphhar, den_mt, nsp, n, grad, ch) + USE m_grdchlh + USE m_mkgylm + IMPLICIT NONE + CLASS(t_xcpot), INTENT(IN) :: xcpot + TYPE(t_atoms), INTENT(IN) :: atoms + TYPE(t_sphhar), INTENT(IN) :: sphhar + REAL, INTENT(IN) :: den_mt(:, 0:, :) + INTEGER, INTENT(IN) :: n, jspins, nsp + REAL, INTENT(OUT), OPTIONAL :: ch(:, :) + TYPE(t_gradients), INTENT(INOUT):: grad + + REAL, ALLOCATABLE :: chlh(:, :, :), chlhdr(:, :, :), chlhdrr(:, :, :) + REAL, ALLOCATABLE :: chdr(:, :), chdt(:, :), chdf(:, :), ch_tmp(:, :) + REAL, ALLOCATABLE :: chdrr(:, :), chdtt(:, :), chdff(:, :), chdtf(:, :) + REAL, ALLOCATABLE :: chdrt(:, :), chdrf(:, :) + INTEGER:: nd, lh, js, jr, kt, k + + nd = atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1) + + ALLOCATE (chlh(atoms%jmtd, 0:sphhar%nlhd, jspins)) + ALLOCATE (ch_tmp(nsp, jspins)) + IF (xcpot%needs_grad()) THEN + ALLOCATE (chdr(nsp, jspins), chdt(nsp, jspins), chdf(nsp, jspins), chdrr(nsp, jspins), & + chdtt(nsp, jspins), chdff(nsp, jspins), chdtf(nsp, jspins), chdrt(nsp, jspins), & + chdrf(nsp, jspins)) + ALLOCATE (chlhdr(atoms%jmtd, 0:sphhar%nlhd, jspins)) + ALLOCATE (chlhdrr(atoms%jmtd, 0:sphhar%nlhd, jspins)) + ENDIF + + DO lh = 0, sphhar%nlh(nd) + + ! calculates gradients of radial charge densities of l=> 0. + ! rho*ylh/r**2 is charge density. chlh=rho/r**2. + ! charge density=sum(chlh*ylh). + ! chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr. + + DO js = 1, jspins + DO jr = 1, atoms%jri(n) + chlh(jr, lh, js) = den_mt(jr, lh, js)/(atoms%rmsh(jr, n)*atoms%rmsh(jr, n)) + ENDDO + IF (xcpot%needs_grad()) CALL grdchlh(1, 1, atoms%jri(n), atoms%dx(n), atoms%rmsh(1, n), & + chlh(1, lh, js), ndvgrd, chlhdr(1, lh, js), chlhdrr(1, lh, js)) + + ENDDO ! js + ENDDO ! lh + + kt = 0 + DO jr = 1, atoms%jri(n) + ! charge density (on extended grid for all jr) + ! following are at points on jr-th sphere. + ch_tmp(:, :) = 0.0 + ! generate the densities on an angular mesh + DO js = 1, jspins + DO lh = 0, sphhar%nlh(nd) + DO k = 1, nsp + ch_tmp(k, js) = ch_tmp(k, js) + ylh(k, lh, nd)*chlh(jr, lh, js) + ENDDO + ENDDO + ENDDO + IF (xcpot%needs_grad()) THEN + chdr(:, :) = 0.0 ! d(ch)/dr + chdt(:, :) = 0.0 ! d(ch)/dtheta + chdf(:, :) = 0.0 ! d(ch)/dfai + chdrr(:, :) = 0.0 ! dd(ch)/drr + chdtt(:, :) = 0.0 ! dd(ch)/dtt + chdff(:, :) = 0.0 ! dd(ch)/dff + chdtf(:, :) = 0.0 ! dd(ch)/dtf + chdrt(:, :) = 0.0 ! d(d(ch)/dr)dt + chdrf(:, :) = 0.0 ! d(d(ch)/dr)df + ! generate the derivatives on an angular mesh + DO js = 1, jspins + DO lh = 0, sphhar%nlh(nd) + ! + DO k = 1, nsp + chdr(k, js) = chdr(k, js) + ylh(k, lh, nd)*chlhdr(jr, lh, js) + chdrr(k, js) = chdrr(k, js) + ylh(k, lh, nd)*chlhdrr(jr, lh, js) + ENDDO + + DO k = 1, nsp + chdrt(k, js) = chdrt(k, js) + ylht(k, lh, nd)*chlhdr(jr, lh, js) + chdrf(k, js) = chdrf(k, js) + ylhf(k, lh, nd)*chlhdr(jr, lh, js) + chdt(k, js) = chdt(k, js) + ylht(k, lh, nd)*chlh(jr, lh, js) + chdf(k, js) = chdf(k, js) + ylhf(k, lh, nd)*chlh(jr, lh, js) + chdtt(k, js) = chdtt(k, js) + ylhtt(k, lh, nd)*chlh(jr, lh, js) + chdff(k, js) = chdff(k, js) + ylhff(k, lh, nd)*chlh(jr, lh, js) + chdtf(k, js) = chdtf(k, js) + ylhtf(k, lh, nd)*chlh(jr, lh, js) + ENDDO + ENDDO ! lh + ENDDO ! js + + CALL mkgylm(jspins, atoms%rmsh(jr, n), thet, nsp, & + ch_tmp, chdr, chdt, chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt) + ENDIF + !Set charge to minimum value + IF (PRESENT(ch)) THEN + WHERE (ABS(ch_tmp) < d_15) ch_tmp = d_15 + ch(kt + 1:kt + nsp, :) = ch_tmp(:nsp, :) + ENDIF + kt = kt + nsp + END DO + + END SUBROUTINE mt_to_grid + + SUBROUTINE mt_from_grid(atoms, sphhar, nsp, n, jspins, v_in, vr) + IMPLICIT NONE + TYPE(t_atoms), INTENT(IN) :: atoms + TYPE(t_sphhar), INTENT(IN):: sphhar + INTEGER, INTENT(IN) :: nsp, jspins, n + REAL, INTENT(IN) :: v_in(:, :) + REAL, INTENT(INOUT) :: vr(:, 0:, :) + + REAL :: vpot(nsp), vlh + INTEGER :: js, kt, lh, jr, nd + nd = atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1) + DO js = 1, jspins + ! + kt = 0 + DO jr = 1, atoms%jri(n) + vpot = v_in(kt + 1:kt + nsp, js)*wt(:)! multiplicate v_in with the weights of the k-points + + DO lh = 0, sphhar%nlh(nd) + ! + ! ---> determine the corresponding potential number + !c through gauss integration + ! + vlh = dot_PRODUCT(vpot(:), ylh(:nsp, lh, nd)) + vr(jr, lh, js) = vr(jr, lh, js) + vlh + ENDDO ! lh + kt = kt + nsp + ENDDO ! jr + ENDDO + + END SUBROUTINE mt_from_grid + + SUBROUTINE finish_mt_grid() + DEALLOCATE (ylh, wt, rx, thet) + IF (ALLOCATED(ylht)) DEALLOCATE (ylht, ylhtt, ylhf, ylhff, ylhtf) + END SUBROUTINE finish_mt_grid END MODULE m_mt_tofrom_grid diff --git a/xc-pot/libxc_postprocess_gga.f90 b/xc-pot/libxc_postprocess_gga.f90 index e353ebcbcd48daaa6947e662fb864fde2df0b71f..97be0932419e61a0ffcccb9bed0ea503de033fae 100644 --- a/xc-pot/libxc_postprocess_gga.f90 +++ b/xc-pot/libxc_postprocess_gga.f90 @@ -5,119 +5,85 @@ !-------------------------------------------------------------------------------- MODULE m_libxc_postprocess_gga CONTAINS - - - SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad) - USE m_mt_tofrom_grid - USE m_types - IMPLICIT NONE - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_sphhar),INTENT(IN) :: sphhar - INTEGER,INTENT(IN) :: n - REAL,INTENT(INOUT) :: v_xc(:,:) - TYPE(t_gradients),INTENT(IN):: grad + SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad) + USE m_mt_tofrom_grid + USE m_types + IMPLICIT NONE + CLASS(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_atoms),INTENT(IN) :: atoms + TYPE(t_sphhar),INTENT(IN) :: sphhar + INTEGER,INTENT(IN) :: n + REAL,INTENT(INOUT) :: v_xc(:,:) + TYPE(t_gradients),INTENT(IN):: grad - INTEGER :: nsp,n_sigma,i - REAL,ALLOCATABLE:: vsigma(:,:),vsigma_mt(:,:,:) - TYPE(t_gradients)::grad_vsigma - - n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise - nsp=SIZE(v_xc,1) !no of points - ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:sphhar%nlhd,n_sigma)) - vsigma_mt=0.0 - vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix - CALL mt_from_grid(atoms,sphhar,nsp/atoms%jmtd,n,n_sigma,vsigma,vsigma_mt) - DO i=1,atoms%jri(n) - vsigma_mt(i,:,:)=vsigma_mt(i,:,:)*atoms%rmsh(i,n)**2 - ENDDO - ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma)) - CALL mt_to_grid(xcpot,n_sigma,atoms,sphhar,vsigma_mt,nsp/atoms%jmtd,n,grad=grad_vsigma) - - CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc) - END SUBROUTINE libxc_postprocess_gga_mt + INTEGER :: nsp,n_sigma,i + REAL,ALLOCATABLE:: vsigma(:,:),vsigma_mt(:,:,:) + TYPE(t_gradients)::grad_vsigma - SUBROUTINE libxc_postprocess_gga_pw(xcpot,stars,cell,v_xc,grad) - USE m_pw_tofrom_grid - USE m_types - IMPLICIT NONE - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_stars),INTENT(IN) :: stars - TYPE(t_cell),INTENT(IN) :: cell - REAL,INTENT(INOUT) :: v_xc(:,:) - TYPE(t_gradients),INTENT(IN):: grad + n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise + nsp=SIZE(v_xc,1) !no of points + ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:sphhar%nlhd,n_sigma)) + vsigma_mt=0.0 + vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix + CALL mt_from_grid(atoms,sphhar,nsp/atoms%jmtd,n,n_sigma,vsigma,vsigma_mt) + DO i=1,atoms%jri(n) + vsigma_mt(i,:,:)=vsigma_mt(i,:,:)*atoms%rmsh(i,n)**2 + ENDDO + ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma)) + CALL mt_to_grid(xcpot,n_sigma,atoms,sphhar,vsigma_mt,nsp/atoms%jmtd,n,grad=grad_vsigma) - COMPLEX,ALLOCATABLE:: vsigma_g(:,:) - REAL,ALLOCATABLE:: vsigma(:,:) - TYPE(t_gradients)::grad_vsigma - INTEGER :: nsp,n_sigma - - nsp=SIZE(v_xc,1) !no of points - n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine - ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma));vsigma_g=0.0 - vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix - CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g) - !vsigma_g(:,1)=vsigma_g(:,1)*stars%nstr(:) - ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma)) - CALL pw_to_grid(xcpot,n_sigma,.false.,stars,cell,vsigma_g,grad_vsigma) - - CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc) + CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc) + END SUBROUTINE libxc_postprocess_gga_mt - END SUBROUTINE libxc_postprocess_gga_pw + SUBROUTINE libxc_postprocess_gga_pw(xcpot,stars,cell,v_xc,grad) + USE m_pw_tofrom_grid + USE m_types + IMPLICIT NONE + CLASS(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_stars),INTENT(IN) :: stars + TYPE(t_cell),INTENT(IN) :: cell + REAL,INTENT(INOUT) :: v_xc(:,:) + TYPE(t_gradients),INTENT(IN):: grad + COMPLEX,ALLOCATABLE:: vsigma_g(:,:) + REAL,ALLOCATABLE:: vsigma(:,:) + TYPE(t_gradients)::grad_vsigma + INTEGER :: nsp,n_sigma - SUBROUTINE libxc_postprocess_gga_pw_alt(xcpot,stars,cell,v_xc,grad) - USE m_pw_tofrom_grid - USE m_types - IMPLICIT NONE - CLASS(t_xcpot),INTENT(IN) :: xcpot - TYPE(t_stars),INTENT(IN) :: stars - TYPE(t_cell),INTENT(IN) :: cell - REAL,INTENT(INOUT) :: v_xc(:,:) - TYPE(t_gradients),INTENT(IN):: grad + nsp=SIZE(v_xc,1) !no of points + n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine + ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma)); vsigma_g=0.0 + vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix + CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g) + !vsigma_g(:,1)=vsigma_g(:,1)*stars%nstr(:) + ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma)) + CALL pw_to_grid(xcpot,n_sigma,.false.,stars,cell,vsigma_g,grad_vsigma) - COMPLEX,ALLOCATABLE:: vsigma_g(:,:) - REAL,ALLOCATABLE:: vsigma(:,:) - TYPE(t_gradients)::grad_vsigma - INTEGER :: nsp,n_sigma,i - - nsp=SIZE(v_xc,1) !no of points - n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine - ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma));vsigma_g=0.0 - ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma)) - DO i=1,3 - vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix - !Multiply with gradient - vsigma(:,1)=vsigma(:,1)*grad%gr(i,:,1) - CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g) - CALL pw_to_grid(xcpot,n_sigma,.FALSE.,stars,cell,vsigma_g,grad_vsigma) - v_xc(:,1)=v_xc(:,1)-2*grad_vsigma%gr(i,:,1) - ENDDO - - END SUBROUTINE libxc_postprocess_gga_pw_alt - - SUBROUTINE libxc_postprocess_gga(vsigma,grad,grad_vsigma,v_xc) - USE m_types - IMPLICIT NONE - REAL,INTENT(IN) :: vsigma(:,:) - TYPE(t_gradients),INTENT(IN):: grad,grad_vsigma - REAL,INTENT(INOUT) :: v_xc(:,:) - INTEGER:: i - IF (SIZE(v_xc,2)==1) THEN !Single spin - DO i=1,SIZE(v_xc,1) !loop over points - v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1) - ENDDO - ELSE !two spins - DO i=1,SIZE(v_xc,1) !loop over points - v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)-& - dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,2))-vsigma(i,2)*grad%laplace(i,2) - v_xc(i,2)=v_xc(i,2)-2*dot_PRODUCT(grad_vsigma%gr(:,i,3),grad%gr(:,i,2))-2*vsigma(i,3)*grad%laplace(i,2)-& - dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,1))-vsigma(i,2)*grad%laplace(i,1) - ENDDO - END IF - - END SUBROUTINE libxc_postprocess_gga + CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc) + + END SUBROUTINE libxc_postprocess_gga_pw + + SUBROUTINE libxc_postprocess_gga(vsigma,grad,grad_vsigma,v_xc) + USE m_types + IMPLICIT NONE + REAL,INTENT(IN) :: vsigma(:,:) + TYPE(t_gradients),INTENT(IN):: grad,grad_vsigma + REAL,INTENT(INOUT) :: v_xc(:,:) + INTEGER:: i + IF (SIZE(v_xc,2)==1) THEN !Single spin + DO i=1,SIZE(v_xc,1) !loop over points + v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1) + ENDDO + ELSE !two spins + DO i=1,SIZE(v_xc,1) !loop over points + v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)-& + dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,2))-vsigma(i,2)*grad%laplace(i,2) + v_xc(i,2)=v_xc(i,2)-2*dot_PRODUCT(grad_vsigma%gr(:,i,3),grad%gr(:,i,2))-2*vsigma(i,3)*grad%laplace(i,2)-& + dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,1))-vsigma(i,2)*grad%laplace(i,1) + ENDDO + END IF + + END SUBROUTINE libxc_postprocess_gga - END MODULE m_libxc_postprocess_gga