From d87452fcbca9c7fbeae5848fba98b613b8a88c34 Mon Sep 17 00:00:00 2001 From: ceresoli Date: Mon, 19 Mar 2007 09:38:11 +0000 Subject: [PATCH] GIPAW should not depend on PH! I've just borrowed few routines. That was the origin of a lot of troubles. The routines in PH use the PHCOM module, which was vastly *uninitialized*! We should instead use the GIPAW_MODULE module. Please, do not tie us again with PHONON! (D.C.) git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3858 c92efa57-630b-4861-b058-cf58834340f0 --- GIPAW/Makefile | 14 +- GIPAW/cg_psi.f90 | 41 ++++ GIPAW/cgsolve_all.f90 | 238 ++++++++++++++++++++++ GIPAW/{ch_psi_all2.f90 => ch_psi_all.f90} | 66 ++---- GIPAW/h_psiq.f90 | 102 ++++++++++ GIPAW/h_psiq2.f90 | 157 -------------- GIPAW/make.depend | 27 +-- 7 files changed, 422 insertions(+), 223 deletions(-) create mode 100644 GIPAW/cg_psi.f90 create mode 100644 GIPAW/cgsolve_all.f90 rename GIPAW/{ch_psi_all2.f90 => ch_psi_all.f90} (57%) create mode 100644 GIPAW/h_psiq.f90 delete mode 100644 GIPAW/h_psiq2.f90 diff --git a/GIPAW/Makefile b/GIPAW/Makefile index 163715644..889189535 100644 --- a/GIPAW/Makefile +++ b/GIPAW/Makefile @@ -9,10 +9,12 @@ stop_code.o \ apply_p.o \ apply_vel.o \ greenfunction.o \ -h_psiq2.o \ +h_psiq.o \ +cg_psi.o \ +cgsolve_all.o \ sym_cart_tensor.o \ symmetrize_field.o \ -ch_psi_all2.o \ +ch_psi_all.o \ test_sum_rule.o \ compute_u_kq.o \ suscept_crystal.o \ @@ -89,15 +91,13 @@ MODULES = \ PWOBJS = ../PW/libpw.a -PHOBJS = ../PH/libph.a - -TLDEPS=bindir ph pw mods libs libiotk +TLDEPS=bindir pw mods libs libiotk all : tldeps gipaw.x -gipaw.x : $(GIPAWOBJS) $(PWOBJS) $(PHOBJS) $(LIBOBJS) +gipaw.x : $(GIPAWOBJS) $(PWOBJS) $(LIBOBJS) $(MPIF90) $(LDFLAGS) -o $@ \ - $(GIPAWOBJS) $(MODULES) $(PHOBJS) $(PWOBJS) $(LIBOBJS) $(LIBS) + $(GIPAWOBJS) $(MODULES) $(PWOBJS) $(LIBOBJS) $(LIBS) - ( cd ../bin; ln -fs ../GIPAW/$@ . ) tldeps: diff --git a/GIPAW/cg_psi.f90 b/GIPAW/cg_psi.f90 new file mode 100644 index 000000000..c0011bb38 --- /dev/null +++ b/GIPAW/cg_psi.f90 @@ -0,0 +1,41 @@ +! +! Copyright (C) 2001 PWSCF group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +!----------------------------------------------------------------- +subroutine cg_psi (lda, n, m, psi, h_diag) + !----------------------------------------------------------------- + ! + ! This routine gives a preconditioning to the linear system solver. + ! The preconditioning is diagonal in reciprocal space + ! + ! + USE kinds, only : DP + implicit none + + integer :: lda, n, m + ! input: the leading dimension of the psi vector + ! input: the real dimension of the vector + ! input: the number of vectors + + complex(DP) :: psi (lda, m) + ! inp/out: the vector to be preconditioned + + real(DP) :: h_diag (lda, m) + ! input: the preconditioning vector + + integer :: k, i + ! counter on bands + ! counter on the elements of the vector + ! + do k = 1, m + do i = 1, n + psi (i, k) = psi (i, k) * h_diag (i, k) + enddo + enddo + return +end subroutine cg_psi diff --git a/GIPAW/cgsolve_all.f90 b/GIPAW/cgsolve_all.f90 new file mode 100644 index 000000000..6dc81d683 --- /dev/null +++ b/GIPAW/cgsolve_all.f90 @@ -0,0 +1,238 @@ +! +! Copyright (C) 2001-2004 PWSCF group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +#include "f_defs.h" +!---------------------------------------------------------------------- +subroutine cgsolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, & + ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd) + !---------------------------------------------------------------------- + ! + ! iterative solution of the linear system: + ! + ! ( h - e + Q ) * dpsi = d0psi (1) + ! + ! where h is a complex hermitean matrix, e is a real sca + ! dpsi and d0psi are complex vectors + ! + ! on input: + ! h_psi EXTERNAL name of a subroutine: + ! h_psi(ndim,psi,psip) + ! Calculates H*psi products. + ! Vectors psi and psip should be dimensined + ! (ndmx,nvec). nvec=1 is used! + ! + ! cg_psi EXTERNAL name of a subroutine: + ! g_psi(ndmx,ndim,notcnv,psi,e) + ! which calculates (h-e)^-1 * psi, with + ! some approximation, e.g. (diag(h)-e) + ! + ! e real unperturbed eigenvalue. + ! + ! dpsi contains an estimate of the solution + ! vector. + ! + ! d0psi contains the right hand side vector + ! of the system. + ! + ! ndmx integer row dimension of dpsi, ecc. + ! + ! ndim integer actual row dimension of dpsi + ! + ! ethr real convergence threshold. solution + ! improvement is stopped when the error in + ! eq (1), defined as l.h.s. - r.h.s., becomes + ! less than ethr in norm. + ! + ! on output: dpsi contains the refined estimate of the + ! solution vector. + ! + ! d0psi is corrupted on exit + ! + ! revised (extensively) 6 Apr 1997 by A. Dal Corso & F. Mauri + ! revised (to reduce memory) 29 May 2004 by S. de Gironcoli + ! +#include "f_defs.h" + USE kinds, only : DP + implicit none + ! + ! first the I/O variables + ! + integer :: ndmx, & ! input: the maximum dimension of the vectors + ndim, & ! input: the actual dimension of the vectors + kter, & ! output: counter on iterations + nbnd, & ! input: the number of bands + ik ! input: the k point + + real(DP) :: & + e(nbnd), & ! input: the actual eigenvalue + anorm, & ! output: the norm of the error in the solution + h_diag(ndmx,nbnd), & ! input: an estimate of ( H - \epsilon ) + ethr ! input: the required precision + + complex(DP) :: & + dpsi (ndmx, nbnd), & ! output: the solution of the linear syst + d0psi (ndmx, nbnd) ! input: the known term + + logical :: conv_root ! output: if true the root is converged + external h_psi, & ! input: the routine computing h_psi + cg_psi ! input: the routine computing cg_psi + ! + ! here the local variables + ! + integer, parameter :: maxter = 200 + ! the maximum number of iterations + integer :: iter, ibnd, lbnd + ! counters on iteration, bands + integer , allocatable :: conv (:) + ! if 1 the root is converged + + complex(DP), allocatable :: g (:,:), t (:,:), h (:,:), hold (:,:) + ! the gradient of psi + ! the preconditioned gradient + ! the delta gradient + ! the conjugate gradient + ! work space + complex(DP) :: dcgamma, dclambda, ZDOTC + ! the ratio between rho + ! step length + ! the scalar product + real(DP), allocatable :: rho (:), rhoold (:), eu (:), a(:), c(:) + ! the residue + ! auxiliary for h_diag + real(DP) :: kter_eff + ! account the number of iterations with b + ! coefficient of quadratic form + ! + call start_clock ('cgsolve') + allocate ( g(ndmx,nbnd), t(ndmx,nbnd), h(ndmx,nbnd), hold(ndmx ,nbnd) ) + allocate (a(nbnd), c(nbnd)) + allocate (conv ( nbnd)) + allocate (rho(nbnd),rhoold(nbnd)) + allocate (eu ( nbnd)) + ! WRITE( stdout,*) g,t,h,hold + + kter_eff = 0.d0 + do ibnd = 1, nbnd + conv (ibnd) = 0 + enddo + do iter = 1, maxter + ! + ! compute the gradient. can reuse information from previous step + ! + if (iter == 1) then + call h_psi (ndim, dpsi, g, e, ik, nbnd) + do ibnd = 1, nbnd + call ZAXPY (ndim, (-1.d0,0.d0), d0psi(1,ibnd), 1, g(1,ibnd), 1) + enddo + endif + ! + ! compute preconditioned residual vector and convergence check + ! + lbnd = 0 + do ibnd = 1, nbnd + if (conv (ibnd) .eq.0) then + lbnd = lbnd+1 + call ZCOPY (ndim, g (1, ibnd), 1, h (1, ibnd), 1) + call cg_psi(ndmx, ndim, 1, h(1,ibnd), h_diag(1,ibnd) ) + rho(lbnd) = ZDOTC (ndim, h(1,ibnd), 1, g(1,ibnd), 1) + endif + enddo + kter_eff = kter_eff + DBLE (lbnd) / DBLE (nbnd) +#ifdef __PARA + call reduce (lbnd, rho ) +#endif + do ibnd = nbnd, 1, -1 + if (conv(ibnd).eq.0) then + rho(ibnd)=rho(lbnd) + lbnd = lbnd -1 + anorm = sqrt (rho (ibnd) ) + if (anorm.lt.ethr) conv (ibnd) = 1 + endif + enddo +! + conv_root = .true. + do ibnd = 1, nbnd + conv_root = conv_root.and. (conv (ibnd) .eq.1) + enddo + if (conv_root) goto 100 + ! + ! compute the step direction h. Conjugate it to previous step + ! + lbnd = 0 + do ibnd = 1, nbnd + if (conv (ibnd) .eq.0) then +! +! change sign to h +! + call DSCAL (2 * ndim, - 1.d0, h (1, ibnd), 1) + if (iter.ne.1) then + dcgamma = rho (ibnd) / rhoold (ibnd) + call ZAXPY (ndim, dcgamma, hold (1, ibnd), 1, h (1, ibnd), 1) + endif + +! +! here hold is used as auxiliary vector in order to efficiently compute t = A*h +! it is later set to the current (becoming old) value of h +! + lbnd = lbnd+1 + call ZCOPY (ndim, h (1, ibnd), 1, hold (1, lbnd), 1) + eu (lbnd) = e (ibnd) + endif + enddo + ! + ! compute t = A*h + ! + call h_psi (ndim, hold, t, eu, ik, lbnd) + ! + ! compute the coefficients a and c for the line minimization + ! compute step length lambda + lbnd=0 + do ibnd = 1, nbnd + if (conv (ibnd) .eq.0) then + lbnd=lbnd+1 + a(lbnd) = ZDOTC (ndim, h(1,ibnd), 1, g(1,ibnd), 1) + c(lbnd) = ZDOTC (ndim, h(1,ibnd), 1, t(1,lbnd), 1) + end if + end do +#ifdef __PARA + call reduce (lbnd, a) + call reduce (lbnd, c) +#endif + lbnd=0 + do ibnd = 1, nbnd + if (conv (ibnd) .eq.0) then + lbnd=lbnd+1 + dclambda = CMPLX ( - a(lbnd) / c(lbnd), 0.d0) + ! + ! move to new position + ! + call ZAXPY (ndim, dclambda, h(1,ibnd), 1, dpsi(1,ibnd), 1) + ! + ! update to get the gradient + ! + !g=g+lam + call ZAXPY (ndim, dclambda, t(1,lbnd), 1, g(1,ibnd), 1) + ! + ! save current (now old) h and rho for later use + ! + call ZCOPY (ndim, h(1,ibnd), 1, hold(1,ibnd), 1) + rhoold (ibnd) = rho (ibnd) + endif + enddo + enddo +100 continue + kter = kter_eff + deallocate (eu) + deallocate (rho, rhoold) + deallocate (conv) + deallocate (a,c) + deallocate (g, t, h, hold) + + call stop_clock ('cgsolve') + return +end subroutine cgsolve_all diff --git a/GIPAW/ch_psi_all2.f90 b/GIPAW/ch_psi_all.f90 similarity index 57% rename from GIPAW/ch_psi_all2.f90 rename to GIPAW/ch_psi_all.f90 index c8534328a..6b90eb2d1 100644 --- a/GIPAW/ch_psi_all2.f90 +++ b/GIPAW/ch_psi_all.f90 @@ -7,7 +7,7 @@ ! !----------------------------------------------------------------------- -subroutine ch_psi_all2 (n, h, ah, e, ik, m) +subroutine ch_psi_all (n, h, ah, e, ik, m) !----------------------------------------------------------------------- ! ! This routine applies the operator ( H - \epsilon S + alpha_pv P_v) @@ -15,11 +15,10 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m) ! #include "f_defs.h" - use pwcom - use becmod - USE noncollin_module, ONLY : noncolin, npol + USE pwcom + USE becmod USE kinds, only : DP - use phcom + USE gipaw_module implicit none integer :: n, m, ik @@ -30,7 +29,7 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m) real(DP) :: e (m) ! input: the eigenvalue - complex(DP) :: h (npwx*npol, m), ah (npwx*npol, m) + complex(DP) :: h (npwx, m), ah (npwx, m) ! input: the vector ! output: the operator applied to the vector ! @@ -48,81 +47,56 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m) call start_clock ('ch_psi') allocate (ps ( nbnd , m)) - allocate (hpsi( npwx*npol , m)) - allocate (spsi( npwx*npol , m)) + allocate (hpsi( npwx , m)) + allocate (spsi( npwx , m)) hpsi (:,:) = (0.d0, 0.d0) spsi (:,:) = (0.d0, 0.d0) ! ! compute the product of the hamiltonian with the h vector ! - call h_psiq2 (npwx, n, m, h, hpsi, spsi) + call h_psiq (npwx, n, m, h, hpsi, spsi) call start_clock ('last') ! ! then we compute the operator H-epsilon S ! - ah=(0.d0,0.d0) do ibnd = 1, m do ig = 1, n ah (ig, ibnd) = hpsi (ig, ibnd) - e (ibnd) * spsi (ig, ibnd) enddo enddo - IF (noncolin) THEN - do ibnd = 1, m - do ig = 1, n - ah (ig+npwx,ibnd)=hpsi(ig+npwx,ibnd)-e(ibnd)*spsi(ig+npwx,ibnd) - end do - end do - END IF ! ! Here we compute the projector in the valence band ! - ikq = ik + !!!if (lgamma) then + ikq = ik + !!!else + !!! ikq = 2 * ik + !!!endif ps (:,:) = (0.d0, 0.d0) - IF (noncolin) THEN - call ZGEMM ('C', 'N', nbnd_occ (ikq) , m, npwx*npol, (1.d0, 0.d0) , evq, & - npwx*npol, spsi, npwx*npol, (0.d0, 0.d0) , ps, nbnd) - ELSE - call ZGEMM ('C', 'N', nbnd_occ (ikq) , m, n, (1.d0, 0.d0) , evq, & + call ZGEMM ('C', 'N', nbnd_occ (ikq) , m, n, (1.d0, 0.d0) , evq, & npwx, spsi, npwx, (0.d0, 0.d0) , ps, nbnd) - ENDIF ps (:,:) = ps(:,:) * alpha_pv #ifdef __PARA call reduce (2 * nbnd * m, ps) #endif hpsi (:,:) = (0.d0, 0.d0) - IF (noncolin) THEN - call ZGEMM ('N', 'N', npwx*npol, m, nbnd_occ (ikq) , (1.d0, 0.d0) , evq, & - npwx*npol, ps, nbnd, (1.d0, 0.d0) , hpsi, npwx*npol) - ELSE - call ZGEMM ('N', 'N', n, m, nbnd_occ (ikq) , (1.d0, 0.d0) , evq, & - npwx, ps, nbnd, (1.d0, 0.d0) , hpsi, npwx) - END IF + call ZGEMM ('N', 'N', n, m, nbnd_occ (ikq) , (1.d0, 0.d0) , evq, & + npwx, ps, nbnd, (1.d0, 0.d0) , hpsi, npwx) spsi(:,:) = hpsi(:,:) ! ! And apply S again ! - IF (noncolin) THEN - call ccalbec_nc (nkb, npwx, n, npol, m, becp_nc, vkb, hpsi) - call s_psi_nc (npwx, n, m, hpsi, spsi) - ELSE - call ccalbec (nkb, npwx, n, m, becp, vkb, hpsi) - call s_psi (npwx, n, m, hpsi, spsi) - END IF + call ccalbec (nkb, npwx, n, m, becp, vkb, hpsi) + + call s_psi (npwx, n, m, hpsi, spsi) do ibnd = 1, m do ig = 1, n ah (ig, ibnd) = ah (ig, ibnd) + spsi (ig, ibnd) enddo enddo - IF (noncolin) THEN - do ibnd = 1, m - do ig = 1, n - ah (ig+npwx, ibnd) = ah (ig+npwx, ibnd) + spsi (ig+npwx, ibnd) - enddo - enddo - END IF deallocate (spsi) deallocate (hpsi) @@ -130,4 +104,4 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m) call stop_clock ('last') call stop_clock ('ch_psi') return -end subroutine ch_psi_all2 +end subroutine ch_psi_all diff --git a/GIPAW/h_psiq.f90 b/GIPAW/h_psiq.f90 new file mode 100644 index 000000000..79a199a33 --- /dev/null +++ b/GIPAW/h_psiq.f90 @@ -0,0 +1,102 @@ +! +! Copyright (C) 2001 PWSCF group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +!----------------------------------------------------------------------- +subroutine h_psiq (lda, n, m, psi, hpsi, spsi) + !----------------------------------------------------------------------- + ! + ! This routine computes the product of the Hamiltonian + ! and of the S matrix with a m wavefunctions contained + ! in psi. It first computes the bec matrix of these + ! wavefunctions and then with the routines hus_1psi and + ! s_psi computes for each band the required products + ! + + use pwcom + USE wavefunctions_module, ONLY: psic + USE becmod, ONLY: becp + USE kinds, only : DP + use gipaw_module + implicit none + ! + ! Here the local variables + ! + integer :: ibnd + ! counter on bands + + integer :: lda, n, m + ! input: the leading dimension of the array psi + ! input: the real dimension of psi + ! input: the number of psi to compute + integer :: j + ! do loop index + + complex(DP) :: psi (lda, m), hpsi (lda, m), spsi (lda, m) + ! input: the functions where to apply H and S + ! output: H times psi + ! output: S times psi (Us PP's only) + + + call start_clock ('h_psiq') + call start_clock ('init') + + call ccalbec (nkb, npwx, n, m, becp, vkb, psi) + ! + ! Here we apply the kinetic energy (k+G)^2 psi + ! + do ibnd = 1, m + do j = 1, n + hpsi (j, ibnd) = g2kin (j) * psi (j, ibnd) + enddo + enddo + call stop_clock ('init') + ! + ! the local potential V_Loc psi. First the psi in real space + ! + + do ibnd = 1, m + call start_clock ('firstfft') + psic(:) = (0.d0, 0.d0) + !do j = 1, n + ! psic (nls(igk(j))) = psi (j, ibnd) + !enddo + psic (nls(igk(1:n))) = psi(1:n, ibnd) + call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2) + call stop_clock ('firstfft') + ! + ! and then the product with the potential vrs = (vltot+vr) on the smoo + ! + !do j = 1, nrxxs + ! psic (j) = psic (j) * vrs (j, current_spin) + !enddo + psic (1:nrxxs) = psic (1:nrxxs) * vrs (1:nrxxs, current_spin) + ! + ! back to reciprocal space + ! + call start_clock ('secondfft') + call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2) + ! + ! addition to the total product + ! + !do j = 1, n + ! hpsi (j, ibnd) = hpsi (j, ibnd) + psic (nls(igk(j))) + !enddo + hpsi (1:n, ibnd) = hpsi (1:n, ibnd) + psic (nls(igk(1:n))) + call stop_clock ('secondfft') + enddo + ! + ! Here the product with the non local potential V_NL psi + ! + + call add_vuspsi (lda, n, m, psi, hpsi) + + call s_psi (lda, n, m, psi, spsi) + + call stop_clock ('h_psiq') + return +end subroutine h_psiq diff --git a/GIPAW/h_psiq2.f90 b/GIPAW/h_psiq2.f90 deleted file mode 100644 index 40532bde3..000000000 --- a/GIPAW/h_psiq2.f90 +++ /dev/null @@ -1,157 +0,0 @@ -! -! Copyright (C) 2001 PWSCF group -! This file is distributed under the terms of the -! GNU General Public License. See the file `License' -! in the root directory of the present distribution, -! or http://www.gnu.org/copyleft/gpl.txt . -! -! -!----------------------------------------------------------------------- -subroutine h_psiq2 (lda, n, m, psi, hpsi, spsi) - !----------------------------------------------------------------------- - ! - ! This routine computes the product of the Hamiltonian - ! and of the S matrix with a m wavefunctions contained - ! in psi. It first computes the bec matrix of these - ! wavefunctions and then with the routines hus_1psi and - ! s_psi computes for each band the required products - ! - - use pwcom - USE wavefunctions_module, ONLY: psic, psic_nc - USE becmod, ONLY: becp, becp_nc - USE noncollin_module, ONLY : noncolin, npol - USE kinds, only : DP - use phcom - implicit none - ! - ! Here the local variables - ! - integer :: ibnd - ! counter on bands - - integer :: lda, n, m - ! input: the leading dimension of the array psi - ! input: the real dimension of psi - ! input: the number of psi to compute - integer :: j - ! do loop index - - complex(DP) :: psi (lda*npol, m), hpsi (lda*npol, m), spsi (lda*npol, m) - complex(DP) :: sup, sdwn - ! input: the functions where to apply H and S - ! output: H times psi - ! output: S times psi (Us PP's only) - - - call start_clock ('h_psiq') - call start_clock ('init') - - IF (noncolin) THEN - call ccalbec_nc (nkb, npwx, n, npol, m, becp_nc, vkb, psi) - ELSE - call ccalbec (nkb, npwx, n, m, becp, vkb, psi) - END IF - ! - ! Here we apply the kinetic energy (k+G)^2 psi - ! - hpsi=(0.d0,0.d0) - do ibnd = 1, m - do j = 1, n - hpsi (j, ibnd) = g2kin (j) * psi (j, ibnd) - enddo - enddo - IF (noncolin) THEN - DO ibnd = 1, m - DO j = 1, n - hpsi (j+lda, ibnd) = g2kin (j) * psi (j+lda, ibnd) - ENDDO - ENDDO - ENDIF - call stop_clock ('init') - ! - ! the local potential V_Loc psi. First the psi in real space - ! - - do ibnd = 1, m - call start_clock ('firstfft') - IF (noncolin) THEN - psic_nc = (0.d0, 0.d0) - do j = 1, n - psic_nc(nls(igk(j)),1) = psi (j, ibnd) - psic_nc(nls(igk(j)),2) = psi (j+lda, ibnd) - enddo - call cft3s (psic_nc(1,1), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2) - call cft3s (psic_nc(1,2), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2) - ELSE - psic(:) = (0.d0, 0.d0) - do j = 1, n - psic (nls(igk(j))) = psi (j, ibnd) - enddo - call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2) - END IF - call stop_clock ('firstfft') - ! - ! and then the product with the potential vrs = (vltot+vr) on the smoo - ! - if (noncolin) then - if (domag) then - do j=1, nrxxs - sup = psic_nc(j,1) * (vrs(j,1)+vrs(j,4)) + & - psic_nc(j,2) * (vrs(j,2)-(0.d0,1.d0)*vrs(j,3)) - sdwn = psic_nc(j,2) * (vrs(j,1)-vrs(j,4)) + & - psic_nc(j,1) * (vrs(j,2)+(0.d0,1.d0)*vrs(j,3)) - psic_nc(j,1)=sup - psic_nc(j,2)=sdwn - end do - else - do j=1, nrxxs - psic_nc(j,1)=psic_nc(j,1) * vrs(j,1) - psic_nc(j,2)=psic_nc(j,2) * vrs(j,1) - enddo - endif - else - do j = 1, nrxxs - psic (j) = psic (j) * vrs (j, current_spin) - enddo - endif - ! - ! back to reciprocal space - ! - call start_clock ('secondfft') - IF (noncolin) THEN - call cft3s(psic_nc(1,1),nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,-2) - call cft3s(psic_nc(1,2),nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,-2) - ! - ! addition to the total product - ! - do j = 1, n - hpsi (j, ibnd) = hpsi (j, ibnd) + psic_nc (nls(igk(j)), 1) - hpsi (j+lda, ibnd) = hpsi (j+lda, ibnd) + psic_nc (nls(igk(j)), 2) - enddo - ELSE - call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2) - ! - ! addition to the total product - ! - do j = 1, n - hpsi (j, ibnd) = hpsi (j, ibnd) + psic (nls(igk(j))) - enddo - END IF - call stop_clock ('secondfft') - enddo - ! - ! Here the product with the non local potential V_NL psi - ! - - IF (noncolin) THEN - call add_vuspsi_nc (lda, n, m, psi, hpsi) - call s_psi_nc (lda, n, m, psi, spsi) - ELSE - call add_vuspsi (lda, n, m, psi, hpsi) - call s_psi (lda, n, m, psi, spsi) - END IF - - call stop_clock ('h_psiq') - return -end subroutine h_psiq2 diff --git a/GIPAW/make.depend b/GIPAW/make.depend index cfa61a1ad..0194491ef 100644 --- a/GIPAW/make.depend +++ b/GIPAW/make.depend @@ -10,11 +10,12 @@ biot_savart.o : ../Modules/constants.o biot_savart.o : ../Modules/kind.o biot_savart.o : ../PW/pwcom.o biot_savart.o : gipaw_module.o -ch_psi_all2.o : ../Modules/kind.o -ch_psi_all2.o : ../PH/phcom.o -ch_psi_all2.o : ../PW/becmod.o -ch_psi_all2.o : ../PW/noncol.o -ch_psi_all2.o : ../PW/pwcom.o +cg_psi.o : ../Modules/kind.o +cgsolve_all.o : ../Modules/kind.o +ch_psi_all.o : ../Modules/kind.o +ch_psi_all.o : ../PW/becmod.o +ch_psi_all.o : ../PW/pwcom.o +ch_psi_all.o : gipaw_module.o compute_sigma.o : ../Modules/io_global.o compute_sigma.o : ../Modules/ions_base.o compute_sigma.o : ../Modules/kind.o @@ -75,9 +76,9 @@ gipaw_module.o : ../Modules/mp.o gipaw_module.o : ../Modules/parameters.o gipaw_module.o : ../Modules/splinelib.o gipaw_module.o : ../Modules/uspp.o -gipaw_module.o : ../PW/becmod.o gipaw_module.o : ../PW/paw.o gipaw_module.o : ../PW/pwcom.o +greenfunction.o : ../Modules/io_files.o greenfunction.o : ../Modules/io_global.o greenfunction.o : ../Modules/kind.o greenfunction.o : ../Modules/wavefunctions.o @@ -85,12 +86,11 @@ greenfunction.o : ../PW/becmod.o greenfunction.o : ../PW/noncol.o greenfunction.o : ../PW/pwcom.o greenfunction.o : gipaw_module.o -h_psiq2.o : ../Modules/kind.o -h_psiq2.o : ../Modules/wavefunctions.o -h_psiq2.o : ../PH/phcom.o -h_psiq2.o : ../PW/becmod.o -h_psiq2.o : ../PW/noncol.o -h_psiq2.o : ../PW/pwcom.o +h_psiq.o : ../Modules/kind.o +h_psiq.o : ../Modules/wavefunctions.o +h_psiq.o : ../PW/becmod.o +h_psiq.o : ../PW/pwcom.o +h_psiq.o : gipaw_module.o init_paw_2_no_phase.o : ../Modules/cell_base.o init_paw_2_no_phase.o : ../Modules/constants.o init_paw_2_no_phase.o : ../Modules/ions_base.o @@ -147,7 +147,8 @@ write_tensor_field.o : ../Modules/mp_global.o write_tensor_field.o : ../PW/pwcom.o write_tensor_field.o : gipaw_module.o apply_vel.o : ../include/f_defs.h -ch_psi_all2.o : ../include/f_defs.h +cgsolve_all.o : ../include/f_defs.h +ch_psi_all.o : ../include/f_defs.h compute_u_kq.o : ../include/f_defs.h efg.o : ../include/f_defs.h g_tensor_crystal.o : ../include/f_defs.h