Commit e99aaa5f authored by Gregor Michalicek's avatar Gregor Michalicek

Introduce t_force type to cdn/cdnval.F90

parent bf160b72
......@@ -164,16 +164,12 @@ CONTAINS
REAL, ALLOCATABLE :: svac(:,:),pvac(:,:),mcd(:,:,:)
REAL, ALLOCATABLE :: enerlo(:,:,:),qmat(:,:,:,:)
COMPLEX, ALLOCATABLE :: acof(:,:,:,:),bcof(:,:,:,:),ccof(:,:,:,:,:)
COMPLEX, ALLOCATABLE :: acoflo(:,:,:,:),bcoflo(:,:,:,:)
COMPLEX, ALLOCATABLE :: cveccof(:,:,:,:,:),f_a12(:,:)
COMPLEX, ALLOCATABLE :: e1cof(:,:,:),e2cof(:,:,:),f_a21(:,:)
COMPLEX, ALLOCATABLE :: f_b4(:,:),f_b8(:,:)
COMPLEX, ALLOCATABLE :: aveccof(:,:,:,:),bveccof(:,:,:,:)
COMPLEX, ALLOCATABLE :: qstars(:,:,:,:),m_mcd(:,:,:,:)
TYPE (t_orb) :: orb
TYPE (t_denCoeffs) :: denCoeffs
TYPE (t_denCoeffsOffdiag) :: denCoeffsOffdiag
TYPE (t_force) :: force
TYPE (t_usdus) :: usdus
TYPE (t_zMat) :: zMat
......@@ -220,23 +216,15 @@ CONTAINS
CALL usdus%init(atoms,input%jspins)
CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
CALL denCoeffsOffdiag%init(atoms,noco,sphhar,l_fmpl)
CALL force%init1(input,atoms)
IF ((l_fmpl).AND.(.not.noco%l_mperp)) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval")
svac(:,:) = 0.0 ; pvac(:,:) = 0.0
sqal(:,:,:) = 0.0 ; ener(:,:,:) = 0.0
!+soc
CALL orb%init(atoms,noco,jsp_start,jsp_end)
!+for
IF (input%l_f) THEN
ALLOCATE ( f_a12(3,atoms%ntype),f_a21(3,atoms%ntype) )
ALLOCATE ( f_b4(3,atoms%ntype),f_b8(3,atoms%ntype) )
f_b4(:,:) = czero ; f_a12(:,:) = czero
f_b8(:,:) = czero ; f_a21(:,:) = czero
ELSE
ALLOCATE ( f_b8(1,1) )
ENDIF
!
INQUIRE (file='mcd_inp',exist=l_mcd)
IF (l_mcd) THEN
OPEN (23,file='mcd_inp',STATUS='old',FORM='formatted')
......@@ -342,9 +330,8 @@ CONTAINS
ALLOCATE ( nmtsl(atoms%ntype,nsld),nslat(atoms%nat,nsld) )
ALLOCATE ( zsl(2,nsld),volsl(nsld) )
ALLOCATE ( volintsl(nsld) )
CALL slabgeom(&
atoms,cell,nsld,&
nsl,zsl,nmtsl,nslat,volsl,volintsl)
CALL slabgeom(atoms,cell,nsld,&
nsl,zsl,nmtsl,nslat,volsl,volintsl)
ALLOCATE ( qintsl(nsld,dimension%neigd))
ALLOCATE ( qmtsl(nsld,dimension%neigd))
......@@ -556,7 +543,7 @@ CONTAINS
IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
CALL timestart("cdnval: pwden")
CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
jspin,lapw,noccbd,we,eig,den,qis,results%force,f_b8,zMat)
jspin,lapw,noccbd,we,eig,den,qis,results,force%f_b8,zMat)
CALL timestop("cdnval: pwden")
END IF
!+new
......@@ -605,17 +592,10 @@ CONTAINS
DO ispin = jsp_start,jsp_end
IF (input%l_f) THEN
CALL timestart("cdnval: to_pulay")
ALLOCATE (e1cof(noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),&
! Deallocated after call to force_a21
e2cof(noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),&
acoflo(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat),&
bcoflo(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat),&
aveccof(3,noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),&
bveccof(3,noccbd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat),&
cveccof(3,-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat) )
CALL force%init2(noccbd,input,atoms)
CALL abcof(input,atoms,sym, cell,lapw,noccbd,usdus, noco,ispin,oneD,&
acof(:,0:,:,ispin),bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin),zMat,&
eig,acoflo,bcoflo,e1cof,e2cof,aveccof,bveccof,cveccof)
eig,force)
CALL timestop("cdnval: to_pulay")
ELSE
CALL timestart("cdnval: abcof")
......@@ -697,15 +677,11 @@ CONTAINS
IF (.not.input%l_useapw) THEN
CALL force_a12(atoms,noccbd,sym, dimension,cell,oneD,&
we,ispin,noccbd,usdus,acof(:,0:,:,ispin),&
bcof(:,0:,:,ispin),e1cof,e2cof, acoflo,bcoflo, results,f_a12)
bcof(:,0:,:,ispin),force,results)
ENDIF
CALL force_a21(input,atoms,dimension,noccbd,sym,&
oneD,cell,we,ispin,enpara%el0(0:,:,ispin),noccbd,eig,usdus,acof(:,0:,:,ispin),&
bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin), aveccof,bveccof,cveccof,&
results,f_a21,f_b4)
DEALLOCATE (e1cof,e2cof,aveccof,bveccof)
DEALLOCATE (acoflo,bcoflo,cveccof)
bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin),force,results)
CALL timestop("cdnval: force_a12/21")
END IF
......@@ -790,13 +766,11 @@ CONTAINS
CALL timestart("cdnval: dos")
IF (mpi%irank==0) THEN
CALL doswrite(eig_id,dimension,kpts,atoms,vacuum,input,banddos,&
sliceplot,noco,sym,cell,&
l_mcd,ncored,ncore,e_mcd,&
sliceplot,noco,sym,cell,l_mcd,ncored,ncore,e_mcd,&
results%ef,results%bandgap,nsld,oneD)
IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
CALL Ek_write_sl(eig_id,dimension,kpts,atoms,vacuum,&
nsld,input,jspin,sym,cell,&
nsl,nslat)
nsld,input,jspin,sym,cell,nsl,nslat)
END IF
END IF
#ifdef CPP_MPI
......@@ -809,17 +783,16 @@ CONTAINS
CALL cdnmt(dimension%jspd,atoms,sphhar,llpd,&
noco,l_fmpl,jsp_start,jsp_end,&
enpara%el0,enpara%ello0,vTot%mt(:,0,:,:),denCoeffs,&
usdus,orb,&
denCoeffsOffdiag,&
usdus,orb,denCoeffsOffdiag,&
chmom,clmom,qa21,den%mt)
DO ispin = jsp_start,jsp_end
IF (.NOT.sliceplot%slice) THEN
DO n=1,atoms%ntype
enpara%el1(0:3,n,ispin)=ener(0:3,n,ispin)/sqal(0:3,n,ispin)
IF (atoms%nlo(n)>0) enpara%ello1(:atoms%nlo(n),n,ispin)=enerlo(:atoms%nlo(n),n,ispin)/sqlo(:atoms%nlo(n),n,ispin)
ENDDO
IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=pvac(:vacuum%nvac,ispin)/svac(:vacuum%nvac,ispin)
enpara%el1(0:3,n,ispin)=ener(0:3,n,ispin)/sqal(0:3,n,ispin)
IF (atoms%nlo(n)>0) enpara%ello1(:atoms%nlo(n),n,ispin)=enerlo(:atoms%nlo(n),n,ispin)/sqlo(:atoms%nlo(n),n,ispin)
END DO
IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=pvac(:vacuum%nvac,ispin)/svac(:vacuum%nvac,ispin)
END IF
!---> check continuity of charge density
......@@ -835,8 +808,7 @@ CONTAINS
!---> forces of equ. A8 of Yu et al.
IF ((input%l_f)) THEN
CALL timestart("cdnval: force_a8")
CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,&
f_a12,f_a21,f_b4,f_b8,results%force)
CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
CALL timestop("cdnval: force_a8")
END IF
!-for
......
......@@ -7,7 +7,7 @@
MODULE m_pwden
CONTAINS
SUBROUTINE pwden(stars,kpts,banddos,oneD, input,mpi,noco,cell,atoms,sym, &
ikpt,jspin,lapw,ne,we,eig,den,qis,forces,f_b8,zMat)
ikpt,jspin,lapw,ne,we,eig,den,qis,results,f_b8,zMat)
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! In this subroutine the star function expansion coefficients of
! the plane wave charge density is determined.
......@@ -93,6 +93,7 @@ CONTAINS
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_zMat),INTENT(IN) :: zMat
TYPE(t_potden),INTENT(INOUT) :: den
TYPE(t_results),INTENT(INOUT) :: results
REAL,INTENT(IN) :: we(:) !(nobd)
REAL,INTENT(IN) :: eig(:)!(dimension%neigd)
......@@ -102,7 +103,6 @@ CONTAINS
INTEGER,INTENT(IN) :: ikpt,jspin
REAL,INTENT(OUT) :: qis(:,:,:) !(dimension%neigd,kpts%nkpt,dimension%jspd)
COMPLEX, INTENT (INOUT) :: f_b8(3,atoms%ntype)
REAL, INTENT (INOUT) :: forces(:,:,:) !(3,atoms%ntype,dimension%jspd)
!
!-----> LOCAL VARIABLES
!
......@@ -650,7 +650,7 @@ CONTAINS
DO istr = 1 , stars%ng3_fft
ecwk(istr) = scale * ecwk(istr) / REAL( stars%nstr(istr) )
ENDDO
CALL force_b8(atoms,ecwk,stars, sym,cell, jspin, forces,f_b8)
CALL force_b8(atoms,ecwk,stars, sym,cell, jspin, results%force,f_b8)
ENDIF
ENDIF
!
......
......@@ -2,7 +2,7 @@ MODULE m_abcof
CONTAINS
SUBROUTINE abcof(input,atoms,sym, cell,lapw,ne,usdus,&
noco,jspin,oneD, acof,bcof,ccof,zMat,&
eig,acoflo,bcoflo,e1cof,e2cof,aveccof,bveccof,cveccof)
eig,force)
! ************************************************************
! subroutine constructs the a,b coefficients of the linearized
! m.t. wavefunctions for each band and atom. c.l. fu
......@@ -26,6 +26,7 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_zMat),INTENT(IN) :: zMat
TYPE(t_force),OPTIONAL,INTENT(INOUT) :: force
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ne
......@@ -36,13 +37,6 @@ CONTAINS
COMPLEX, INTENT (OUT) :: bcof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
COMPLEX, INTENT (OUT) :: ccof(-atoms%llod:,:,:,:)!(-llod:llod,nobd,atoms%nlod,atoms%nat)
REAL, OPTIONAL, INTENT (IN) :: eig(:)!(dimension%neigd)
COMPLEX, OPTIONAL, INTENT (OUT) :: acoflo(-atoms%llod:,:,:,:)
COMPLEX, OPTIONAL, INTENT (OUT) :: bcoflo(-atoms%llod:,:,:,:)
COMPLEX, OPTIONAL, INTENT (OUT) :: e1cof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
COMPLEX, OPTIONAL, INTENT (OUT) :: e2cof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
COMPLEX, OPTIONAL, INTENT (OUT) :: aveccof(:,:,0:,:)!(3,nobd,0:dimension%lmd,atoms%nat)
COMPLEX, OPTIONAL, INTENT (OUT) :: bveccof(:,:,0:,:)!(3,nobd,0:dimension%lmd,atoms%nat)
COMPLEX, OPTIONAL, INTENT (OUT) :: cveccof(:,-atoms%llod:,:,:,:)
! ..
! .. Local Scalars ..
COMPLEX cexp,phase,c_0,c_1,c_2,ci
......@@ -75,13 +69,13 @@ CONTAINS
bcof(:,:,:) = CMPLX(0.0,0.0)
ccof(:,:,:,:) = CMPLX(0.0,0.0)
IF(PRESENT(eig)) THEN
acoflo(:,:,:,:) = CMPLX(0.0,0.0)
bcoflo(:,:,:,:) = CMPLX(0.0,0.0)
e1cof(:,:,:) = CMPLX(0.0,0.0)
e2cof(:,:,:) = CMPLX(0.0,0.0)
aveccof(:,:,:,:) = CMPLX(0.0,0.0)
bveccof(:,:,:,:) = CMPLX(0.0,0.0)
cveccof(:,:,:,:,:) = CMPLX(0.0,0.0)
force%acoflo = CMPLX(0.0,0.0)
force%bcoflo = CMPLX(0.0,0.0)
force%e1cof = CMPLX(0.0,0.0)
force%e2cof = CMPLX(0.0,0.0)
force%aveccof = CMPLX(0.0,0.0)
force%bveccof = CMPLX(0.0,0.0)
force%cveccof = CMPLX(0.0,0.0)
END IF
! ..
!+APW_LO
......@@ -121,7 +115,7 @@ CONTAINS
!$OMP& jspin,qss,&
!$OMP& apw,const,&
!$OMP& nbasf0,enough,&
!$OMP& acof,bcof,ccof,e1cof,e2cof,acoflo,bcoflo,aveccof,bveccof,cveccof)
!$OMP& acof,bcof,ccof,force)
DO n = 1,atoms%ntype
CALL setabc1lo(atoms,n,usdus,jspin,alo1,blo1,clo1)
......@@ -246,18 +240,18 @@ CONTAINS
IF ((atoms%l_geo(n)).AND.(PRESENT(eig))) THEN
IF (zmat%l_real) THEN
e1cof(:ne,lm,natom) = e1cof(:ne,lm,natom) + c_1 * work_r(:ne) * s2h_e(:ne)
e2cof(:ne,lm,natom) = e2cof(:ne,lm,natom) + c_2 * work_r(:ne) * s2h_e(:ne)
force%e1cof(:ne,lm,natom) = force%e1cof(:ne,lm,natom) + c_1 * work_r(:ne) * s2h_e(:ne)
force%e2cof(:ne,lm,natom) = force%e2cof(:ne,lm,natom) + c_2 * work_r(:ne) * s2h_e(:ne)
DO i = 1,3
aveccof(i,:ne,lm,natom) = aveccof(i,:ne,lm,natom) + c_1 * work_r(:ne) * fgp(i)
bveccof(i,:ne,lm,natom) = bveccof(i,:ne,lm,natom) + c_2 * work_r(:ne) * fgp(i)
force%aveccof(i,:ne,lm,natom) = force%aveccof(i,:ne,lm,natom) + c_1 * work_r(:ne) * fgp(i)
force%bveccof(i,:ne,lm,natom) = force%bveccof(i,:ne,lm,natom) + c_2 * work_r(:ne) * fgp(i)
END DO
ELSE
e1cof(:ne,lm,natom) = e1cof(:ne,lm,natom) + c_1 * work_c(:ne) * s2h_e(:ne)
e2cof(:ne,lm,natom) = e2cof(:ne,lm,natom) + c_2 * work_c(:ne) * s2h_e(:ne)
force%e1cof(:ne,lm,natom) = force%e1cof(:ne,lm,natom) + c_1 * work_c(:ne) * s2h_e(:ne)
force%e2cof(:ne,lm,natom) = force%e2cof(:ne,lm,natom) + c_2 * work_c(:ne) * s2h_e(:ne)
DO i = 1,3
aveccof(i,:ne,lm,natom) = aveccof(i,:ne,lm,natom) + c_1 * work_c(:ne) * fgp(i)
bveccof(i,:ne,lm,natom) = bveccof(i,:ne,lm,natom) + c_2 * work_c(:ne) * fgp(i)
force%aveccof(i,:ne,lm,natom) = force%aveccof(i,:ne,lm,natom) + c_1 * work_c(:ne) * fgp(i)
force%bveccof(i,:ne,lm,natom) = force%bveccof(i,:ne,lm,natom) + c_2 * work_c(:ne) * fgp(i)
END DO
END IF
ENDIF
......@@ -272,11 +266,11 @@ CONTAINS
CALL CPP_BLAS_caxpy(ne,c_1,work_c,1, acof(1,lmp,jatom),1)
CALL CPP_BLAS_caxpy(ne,c_2,work_c,1, bcof(1,lmp,jatom),1)
IF ((atoms%l_geo(n)).AND.(PRESENT(eig))) THEN
CALL CPP_BLAS_caxpy(ne,c_1,work_c*s2h_e,1, e1cof(1,lmp,jatom),1)
CALL CPP_BLAS_caxpy(ne,c_2,work_c*s2h_e,1, e2cof(1,lmp,jatom),1)
CALL CPP_BLAS_caxpy(ne,c_1,work_c*s2h_e,1, force%e1cof(1,lmp,jatom),1)
CALL CPP_BLAS_caxpy(ne,c_2,work_c*s2h_e,1, force%e2cof(1,lmp,jatom),1)
DO i = 1,3
CALL CPP_BLAS_caxpy(ne,c_1,work_c*fgp(i),1, aveccof(i,1,lmp,jatom),3)
CALL CPP_BLAS_caxpy(ne,c_2,work_c*fgp(i),1, bveccof(i,1,lmp,jatom),3)
CALL CPP_BLAS_caxpy(ne,c_1,work_c*fgp(i),1, force%aveccof(i,1,lmp,jatom),3)
CALL CPP_BLAS_caxpy(ne,c_2,work_c*fgp(i),1, force%bveccof(i,1,lmp,jatom),3)
END DO
END IF
ENDIF
......@@ -288,7 +282,7 @@ CONTAINS
IF (k==lapw%kvec(nkvec,lo,natom)) THEN !check if this k-vector has LO attached
CALL abclocdn(atoms,sym,noco,lapw,cell,ccchi(:,jspin),iintsp,phase,ylm,&
n,natom,k,nkvec,lo,ne,alo1,blo1,clo1,acof,bcof,ccof,zMat,&
fgp,acoflo,bcoflo,aveccof,bveccof,cveccof)
fgp,force%acoflo,force%bcoflo,force%aveccof,force%bveccof,force%cveccof)
ENDIF
ENDDO
END DO
......@@ -334,10 +328,10 @@ CONTAINS
DO ie = 1,ne
ccof(m,ie,ilo,jatom) = inv_f * cexp * CONJG( ccof(-m,ie,ilo,iatom))
IF(PRESENT(eig)) THEN
acoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(acoflo(-m,ie,ilo,iatom))
bcoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(bcoflo(-m,ie,ilo,iatom))
force%acoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(force%acoflo(-m,ie,ilo,iatom))
force%bcoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(force%bcoflo(-m,ie,ilo,iatom))
DO i = 1,3
cveccof(i,m,ie,ilo,jatom) = -inv_f * cexp * CONJG(cveccof(i,-m,ie,ilo,iatom))
force%cveccof(i,m,ie,ilo,jatom) = -inv_f * cexp * CONJG(force%cveccof(i,-m,ie,ilo,iatom))
END DO
END IF
ENDDO
......@@ -355,11 +349,11 @@ CONTAINS
END DO
IF ((atoms%l_geo(n)).AND.(PRESENT(eig))) THEN
DO ie = 1,ne
e1cof(ie,lm,jatom) = inv_f * cexp * CONJG(e1cof(ie,lmp,iatom))
e2cof(ie,lm,jatom) = inv_f * cexp * CONJG(e2cof(ie,lmp,iatom))
force%e1cof(ie,lm,jatom) = inv_f * cexp * CONJG(force%e1cof(ie,lmp,iatom))
force%e2cof(ie,lm,jatom) = inv_f * cexp * CONJG(force%e2cof(ie,lmp,iatom))
DO i = 1,3
aveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(aveccof(i,ie,lmp,iatom))
bveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(bveccof(i,ie,lmp,iatom))
force%aveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(force%aveccof(i,ie,lmp,iatom))
force%bveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(force%bveccof(i,ie,lmp,iatom))
END DO
END DO
END IF
......
......@@ -7,12 +7,12 @@ MODULE m_forcea12
CONTAINS
SUBROUTINE force_a12(&
atoms,nobd,sym, DIMENSION, cell,oneD,&
we,jsp,ne,usdus,acof,bcof,e1cof,e2cof,&
acoflo,bcoflo, results,f_a12)
we,jsp,ne,usdus,acof,bcof,force,results)
USE m_types
USE m_constants
IMPLICIT NONE
TYPE(t_force),INTENT(INOUT) :: force
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
......@@ -29,11 +29,6 @@ CONTAINS
REAL, INTENT (IN) :: we(nobd)
COMPLEX, INTENT (IN) :: acof(nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: bcof(nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: e1cof(nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: e2cof(nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: acoflo(-atoms%llod:atoms%llod,nobd,atoms%nlod,atoms%nat)
COMPLEX, INTENT (IN) :: bcoflo(-atoms%llod:atoms%llod,nobd,atoms%nlod,atoms%nat)
COMPLEX, INTENT (INOUT) :: f_a12(3,atoms%ntype)
! ..
! .. Local Scalars ..
COMPLEX a12,cil1,cil2
......@@ -97,8 +92,8 @@ CONTAINS
DO m1 = -l1,l1
lm1 = l1* (l1+1) + m1
DO ie = 1,ne
acof_flapw(ie,lm1) = acof_flapw(ie,lm1) - acoflo(m1,ie,ilo,natrun)
bcof_flapw(ie,lm1) = bcof_flapw(ie,lm1) - bcoflo(m1,ie,ilo,natrun)
acof_flapw(ie,lm1) = acof_flapw(ie,lm1) - force%acoflo(m1,ie,ilo,natrun)
bcof_flapw(ie,lm1) = bcof_flapw(ie,lm1) - force%bcoflo(m1,ie,ilo,natrun)
ENDDO
ENDDO
ENDDO
......@@ -116,8 +111,8 @@ CONTAINS
DO ie = 1,ne
!
a12 = a12 + CONJG(cil1*&
( acof_flapw(ie,lm1)*usdus%us(l1,n,jsp) + bcof_flapw(ie,lm1)*usdus%uds(l1,n,jsp) ))*cil2*&
( e1cof(ie,lm2,natrun)*usdus%us(l2,n,jsp)+ e2cof(ie,lm2,natrun)*usdus%uds(l2,n,jsp))*we(ie)
(acof_flapw(ie,lm1)*usdus%us(l1,n,jsp) + bcof_flapw(ie,lm1)*usdus%uds(l1,n,jsp) ))*cil2*&
(force%e1cof(ie,lm2,natrun)*usdus%us(l2,n,jsp)+ force%e2cof(ie,lm2,natrun)*usdus%uds(l2,n,jsp))*we(ie)
END DO
aaa(1) = alpha(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1+1)
......@@ -234,7 +229,7 @@ CONTAINS
! is also a solution of Schr. equ. if psi is one.
DO i = 1,3
results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a12(i))
f_a12(i,n) = f_a12(i,n) + forc_a12(i)
force%f_a12(i,n) = force%f_a12(i,n) + forc_a12(i)
END DO
!
! write result moved to force_a8
......
......@@ -3,7 +3,7 @@ CONTAINS
SUBROUTINE force_a21(&
input,atoms,DIMENSION,nobd,sym,oneD,cell,&
we,jsp,epar,ne,eig,usdus,&
acof,bcof,ccof,aveccof,bveccof,cveccof, results,f_a21,f_b4)
acof,bcof,ccof,force,results)
! ************************************************************
! Pulay 2nd and 3rd (A17+A20) term force contribution a la Rici
......@@ -30,6 +30,7 @@ CONTAINS
USE m_constants
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_force),INTENT(INOUT) :: force
TYPE(t_results),INTENT(INOUT):: results
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
......@@ -45,13 +46,9 @@ CONTAINS
! .. Array Arguments ..
REAL, INTENT (IN) :: we(nobd),epar(0:atoms%lmaxd,atoms%ntype)
REAL, INTENT (IN) :: eig(DIMENSION%neigd)
COMPLEX, INTENT (INOUT) :: f_a21(3,atoms%ntype),f_b4(3,atoms%ntype)
COMPLEX, INTENT (IN) :: acof(nobd,0:atoms%lmaxd*(atoms%lmaxd+2) ,atoms%nat)
COMPLEX, INTENT (IN) :: bcof(nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: ccof(-atoms%llod:atoms%llod,nobd,atoms%nlod,atoms%nat)
COMPLEX, INTENT (IN) :: aveccof(3,nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: bveccof(3,nobd,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat )
COMPLEX, INTENT (IN) :: cveccof(3,-atoms%llod:atoms%llod,nobd,atoms%nlod,atoms%nat)
! ..
! .. Local Scalars ..
COMPLEX dtd,dtu,utd,utu
......@@ -150,10 +147,10 @@ CONTAINS
END IF
DO i = 1,3
a21(i,natrun) = a21(i,natrun) + 2.0*&
AIMAG( CONJG(acof(ie,lm1,natrun)) *utu*aveccof(i,ie,lm2,natrun)&
+CONJG(acof(ie,lm1,natrun)) *utd*bveccof(i,ie,lm2,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtu*aveccof(i,ie,lm2,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtd*bveccof(i,ie,lm2,natrun))*we(ie)/atoms%neq(n)
AIMAG( CONJG(acof(ie,lm1,natrun)) *utu*force%aveccof(i,ie,lm2,natrun)&
+CONJG(acof(ie,lm1,natrun)) *utd*force%bveccof(i,ie,lm2,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtu*force%aveccof(i,ie,lm2,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtd*force%bveccof(i,ie,lm2,natrun))*we(ie)/atoms%neq(n)
! END i loop
END DO
END IF
......@@ -180,10 +177,10 @@ CONTAINS
DO i = 1,3
DO natrun = natom,natom + atoms%neq(n) - 1
a21(i,natrun) = a21(i,natrun) + 2.0*&
AIMAG(CONJG(acof(ie,lm1,natrun)) *utu*aveccof(i,ie,lm1,natrun)&
+CONJG(acof(ie,lm1,natrun)) *utd*bveccof(i,ie,lm1,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtu*aveccof(i,ie,lm1,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtd*bveccof(i,ie,lm1,natrun)&
AIMAG(CONJG(acof(ie,lm1,natrun)) *utu*force%aveccof(i,ie,lm1,natrun)&
+CONJG(acof(ie,lm1,natrun)) *utd*force%bveccof(i,ie,lm1,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtu*force%aveccof(i,ie,lm1,natrun)&
+CONJG(bcof(ie,lm1,natrun)) *dtd*force%bveccof(i,ie,lm1,natrun)&
)*we(ie) /atoms%neq(n)
END DO
!
......@@ -200,13 +197,13 @@ CONTAINS
!---> add the local orbital and U contribution to a21
!
CALL force_a21_lo(nobd,atoms,jsp,n,we,eig,ne,&
acof,bcof,ccof,aveccof,bveccof,&
cveccof, tlmplm,usdus, a21)
acof,bcof,ccof,force%aveccof,force%bveccof,&
force%cveccof, tlmplm,usdus, a21)
IF ((atoms%n_u.GT.0).AND.(i_u.LE.atoms%n_u)) THEN
CALL force_a21_U(nobd,atoms,i_u,n,jsp,we,ne,&
usdus,v_mmp,acof,bcof,ccof,&
aveccof,bveccof,cveccof, a21)
force%aveccof,force%bveccof,force%cveccof, a21)
END IF
IF (input%l_useapw) THEN
! -> B4 force
......@@ -221,10 +218,10 @@ CONTAINS
we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
CONJG(acof(ie,lm1,natrun)*usdus%us(l1,n,jsp)&
+bcof(ie,lm1,natrun)*usdus%uds(l1,n,jsp))*&
(aveccof(i,ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
+bveccof(i,ie,lm1,natrun)*usdus%duds(l1,n,jsp) )&
-CONJG(aveccof(i,ie,lm1,natrun)*usdus%us(l1,n,jsp)&
+bveccof(i,ie,lm1,natrun)*usdus%uds(l1,n,jsp) )*&
(force%aveccof(i,ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
+force%bveccof(i,ie,lm1,natrun)*usdus%duds(l1,n,jsp) )&
-CONJG(force%aveccof(i,ie,lm1,natrun)*usdus%us(l1,n,jsp)&
+force%bveccof(i,ie,lm1,natrun)*usdus%uds(l1,n,jsp) )*&
(acof(ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
+bcof(ie,lm1,natrun)*usdus%duds(l1,n,jsp)) )
END DO
......@@ -241,15 +238,15 @@ CONTAINS
we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
CONJG( acof(ie,lm1,natrun)* usdus%us(l1,n,jsp)&
+ bcof(ie,lm1,natrun)* usdus%uds(l1,n,jsp) ) *&
cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp)&
force%cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp)&
+ CONJG(ccof(m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
( aveccof(i,ie,lm1,natrun)* usdus%dus(l1,n,jsp)&
+ bveccof(i,ie,lm1,natrun)* usdus%duds(l1,n,jsp)&
+ cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) ) &
- (CONJG( aveccof(i,ie,lm1,natrun) *usdus%us(l1,n,jsp)&
+ bveccof(i,ie,lm1,natrun) *usdus%uds(l1,n,jsp) ) *&
( force%aveccof(i,ie,lm1,natrun)* usdus%dus(l1,n,jsp)&
+ force%bveccof(i,ie,lm1,natrun)* usdus%duds(l1,n,jsp)&
+ force%cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) ) &
- (CONJG( force%aveccof(i,ie,lm1,natrun) *usdus%us(l1,n,jsp)&
+ force%bveccof(i,ie,lm1,natrun) *usdus%uds(l1,n,jsp) ) *&
ccof(m,ie,lo,natrun) *usdus%dulos(lo,n,jsp)&
+ CONJG(cveccof(i,m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
+ CONJG(force%cveccof(i,m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
( acof(ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
+ bcof(ie,lm1,natrun)*usdus%duds(l1,n,jsp)&
+ ccof(m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) ) ) )
......@@ -356,8 +353,8 @@ CONTAINS
! IS ALSO A SOLUTION OF SCHR. EQU. IF PSI IS ONE.
DO i = 1,3
results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a21(i) + forc_b4(i))
f_a21(i,n) = f_a21(i,n) + forc_a21(i)
f_b4(i,n) = f_b4(i,n) + forc_b4(i)
force%f_a21(i,n) = force%f_a21(i,n) + forc_a21(i)
force%f_b4(i,n) = force%f_b4(i,n) + forc_b4(i)
END DO
!
! write result moved to force_a8
......
......@@ -4,9 +4,7 @@ MODULE m_forcea8
!
! ************************************************************
CONTAINS
SUBROUTINE force_a8(input,atoms,sphhar,&
jsp, vr,rho, f_a12,f_a21,f_b4,f_b8,&
force)
SUBROUTINE force_a8(input,atoms,sphhar,jsp,vr,rho,force,results)
!
USE m_intgr, ONLY : intgr3
USE m_constants,ONLY: pi_const,sfp_const
......@@ -17,16 +15,15 @@ CONTAINS
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_force),INTENT(IN) :: force
TYPE(t_results),INTENT(INOUT) :: results
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp
! ..
! .. Array Arguments ..
COMPLEX, INTENT (IN) :: f_a12(3,atoms%ntype),f_a21(3,atoms%ntype)
COMPLEX, INTENT (IN) :: f_b4(3,atoms%ntype),f_b8(3,atoms%ntype)
REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
REAL, INTENT (IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT (INOUT) :: force(:,:,:) !(3,ntypd,jspd)
! ..
! .. Local Scalars ..
COMPLEX aaa,bbb,ccc,ddd,eee,fff
......@@ -240,7 +237,7 @@ CONTAINS
! sum to existing forces
!
DO i = 1,3
force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a8(i))
results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a8(i))
END DO
!
! write result
......@@ -267,8 +264,8 @@ CONTAINS
IF (atoms%l_geo(n)) THEN
WRITE (6,FMT=8030) n
WRITE (16,FMT=8030) n
WRITE (6,FMT=8040) (f_a12(i,n),i=1,3)
WRITE (16,FMT=8040) (f_a12(i,n),i=1,3)
WRITE (6,FMT=8040) (force%f_a12(i,n),i=1,3)
WRITE (16,FMT=8040) (force%f_a12(i,n),i=1,3)
ENDIF
8030 FORMAT (' FORCES: EQUATION A12 FOR ATOM TYPE',i4)
8040 FORMAT (' FX_A12=',2f10.6,' FY_A12=',2f10.6,' FZ_A12=',2f10.6)
......@@ -280,8 +277,8 @@ CONTAINS
IF (atoms%l_geo(n)) THEN
WRITE (6,FMT=8070) n
WRITE (16,FMT=8070) n
WRITE (6,FMT=8080) (f_b4(i,n),i=1,3)
WRITE (16,FMT=8080) (f_b4(i,n),i=1,3)
WRITE (6,FMT=8080) (force%f_b4(i,n),i=1,3)
WRITE (16,FMT=8080) (force%f_b4(i,n),i=1,3)
ENDIF
8070 FORMAT (' FORCES: EQUATION B4 FOR ATOM TYPE',i4)
8080 FORMAT (' FX_B4=',2f10.6,' FY_B4=',2f10.6,' FZ_B4=',2f10.6)
......@@ -292,8 +289,8 @@ CONTAINS
IF (atoms%l_geo(n)) THEN
WRITE (6,FMT=8090) n
WRITE (16,FMT=8090) n
WRITE (6,FMT=8100) (f_b8(i,n),i=1,3)
WRITE (16,FMT=8100) (f_b8(i,n),i=1,3)
WRITE (6,FMT=8100) (force%f_b8(i,n),i=1,3)
WRITE (16,FMT=8100) (force%f_b8(i,n),i=1,3)
ENDIF
8090 FORMAT (' FORCES: EQUATION B8 FOR ATOM TYPE',i4)
8100 FORMAT (' FX_B8=',2f10.6,' FY_B8=',2f10.6,' FZ_B8=',2f10.6)
......@@ -305,8 +302,8 @@ CONTAINS
IF (atoms%l_geo(n)) THEN
WRITE (6,FMT=8050) n
WRITE (16,FMT=8050) n
WRITE (6,FMT=8060) (f_a21(i,n),i=1,3)
WRITE (16,FMT=8060) (f_a21(i,n),i=1,3)
WRITE (6,FMT=8060) (force%f_a21(i,n),i=1,3)
WRITE (16,FMT=8060) (force%f_a21(i,n),i=1,3)
ENDIF
8050 FORMAT (' FORCES: EQUATION A21 FOR ATOM TYPE',i4)
8060 FORMAT (' FX_A21=',2f10.6,' FY_A21=',2f10.6,' FZ_A21=',2f10.6)
......
......@@ -98,7 +98,27 @@ PRIVATE
PROCEDURE,PASS :: init => denCoeffsOffdiag_init
END TYPE t_denCoeffsOffdiag
PUBLIC t_orb, t_denCoeffs, t_denCoeffsOffdiag
TYPE t_force
COMPLEX, ALLOCATABLE :: f_a12(:,:)
COMPLEX, ALLOCATABLE :: f_a21(:,:)
COMPLEX, ALLOCATABLE :: f_b4(:,:)
COMPLEX, ALLOCATABLE :: f_b8(:,:)
COMPLEX, ALLOCATABLE :: e1cof(:,:,:)
COMPLEX, ALLOCATABLE :: e2cof(:,:,:)
COMPLEX, ALLOCATABLE :: aveccof(:,:,:,:)
COMPLEX, ALLOCATABLE :: bveccof(:,:,:,:)
COMPLEX, ALLOCATABLE :: cveccof(:,:,:,:,:)
COMPLEX, ALLOCATABLE :: acoflo(:,:,:,:)
COMPLEX, ALLOCATABLE :: bcoflo(:,:,:,:)
CONTAINS
PROCEDURE,PASS :: init1 => force_init1
PROCEDURE,PASS :: init2 => force_init2
END TYPE t_force
PUBLIC t_orb, t_denCoeffs, t_denCoeffsOffdiag, t_force
CONTAINS
......@@ -338,4 +358,80 @@ SUBROUTINE denCoeffsOffdiag_init(thisDenCoeffsOffdiag, atoms, noco, sphhar, l_fm
END SUBROUTINE denCoeffsOffdiag_init
SUBROUTINE force_init1(thisForce,input,atoms)
USE m_types_setup
IMPLICIT NONE