From 7f179aeefbcacea1d57d9c52ff4bdd334dc7d149 Mon Sep 17 00:00:00 2001 From: cavazzon Date: Tue, 9 May 2006 11:28:40 +0000 Subject: [PATCH] - Bug fix, computation of empty states with FPMD git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3108 c92efa57-630b-4861-b058-cf58834340f0 --- CPV/emptystates.f90 | 94 +++++++++++++++++++++++++-------------------- CPV/forces.f90 | 51 +++++++++--------------- CPV/ksstates.f90 | 35 ++++++++++------- CPV/optical.f90 | 17 ++++++-- CPV/runcg.f90 | 9 +++-- CPV/runcp.f90 | 19 ++++----- CPV/rundiis.f90 | 17 ++++---- 7 files changed, 132 insertions(+), 110 deletions(-) diff --git a/CPV/emptystates.f90 b/CPV/emptystates.f90 index 526b61d46..f272e70bd 100644 --- a/CPV/emptystates.f90 +++ b/CPV/emptystates.f90 @@ -23,15 +23,15 @@ PRIVATE - INTEGER :: max_emp = 0 ! maximum number of iterations + INTEGER :: max_emp = 0 ! maximum number of iterations REAL(DP) :: ethr_emp ! threshold for convergence REAL(DP) :: delt_emp ! delt for the empty states updating REAL(DP) :: emass_emp ! fictitious mass for the empty states - LOGICAL :: prn_emp = .FALSE. + LOGICAL :: prn_emp = .FALSE. CHARACTER(LEN=256) :: fileempty - LOGICAL :: first = .TRUE. + LOGICAL :: first = .TRUE. INTERFACE EMPTY MODULE PROCEDURE EMPTY_SD @@ -76,14 +76,14 @@ LOGICAL FUNCTION readempty( c_emp, wempt ) -! ... This subroutine reads empty states from unit emptyunit + ! ... This subroutine reads empty states from unit emptyunit - USE wave_types, ONLY: wave_descriptor - USE mp_global, ONLY: me_image, nproc_image, intra_image_comm - USE io_global, ONLY: stdout, ionode, ionode_id - USE mp, ONLY: mp_bcast - USE mp_wave, ONLY: splitwf - USE io_files, ONLY: scradir + USE wave_types, ONLY: wave_descriptor + USE mp_global, ONLY: me_image, nproc_image, intra_image_comm + USE io_global, ONLY: stdout, ionode, ionode_id + USE mp, ONLY: mp_bcast + USE mp_wave, ONLY: splitwf + USE io_files, ONLY: scradir USE reciprocal_vectors, ONLY: ig_l2g IMPLICIT none @@ -97,9 +97,10 @@ INTEGER :: nk, ne(2), ngwm_g, nspin COMPLEX(DP), ALLOCATABLE :: ctmp(:) -! -! ... Subroutine Body -! + + ! + ! ... Subroutine Body + ! IF( wempt%nspin < 1 .OR. wempt%nspin > 2 ) & CALL errore( ' readempty ', ' nspin out of range ', 1 ) @@ -176,7 +177,7 @@ SUBROUTINE writeempty( c_emp, wempt ) -! ... This subroutine writes empty states to unit emptyunit + ! ... This subroutine writes empty states to unit emptyunit USE wave_types, ONLY: wave_descriptor USE mp_global, ONLY: me_image, nproc_image, intra_image_comm @@ -191,9 +192,9 @@ INTEGER :: ig, i, ik, nl, ne(2), ngwm_g, nk, ispin, nspin, ngw LOGICAL :: exst COMPLEX(DP), ALLOCATABLE :: ctmp(:) -! -! ... Subroutine Body -! + ! + ! ... Subroutine Body + ! IF( wempt%nspin < 1 .OR. wempt%nspin > 2 ) & CALL errore( ' writeempty ', ' nspin out of range ', 1 ) @@ -246,33 +247,34 @@ SUBROUTINE gram_empty( ispin, tortho, cf, wfill, ce, wempt ) -! This subroutine orthogonalize the empty states CE to the -! filled states CF using gram-shmitd . -! If TORTHO is FALSE the subroutine orthonormalizes the -! empty states CE and orthogonalize them to the CF states. + ! This subroutine orthogonalize the empty states CE to the + ! filled states CF using gram-shmitd . + ! If TORTHO is FALSE the subroutine orthonormalizes the + ! empty states CE and orthogonalize them to the CF states. USE wave_types, ONLY: wave_descriptor - USE mp, ONLY: mp_sum - USE mp_global, ONLY: nproc_image, intra_image_comm + USE mp, ONLY: mp_sum + USE mp_global, ONLY: nproc_image, intra_image_comm REAL(DP) SQRT, DNRM2 -! ... ARGUMENTS + ! ... ARGUMENTS + LOGICAL, INTENT(IN) :: TORTHO COMPLEX(DP), INTENT(INOUT) :: CF(:,:), CE(:,:) type (wave_descriptor), INTENT(IN) :: wfill, wempt INTEGER, INTENT(IN) :: ispin -! ... LOCALS + ! ... LOCALS + INTEGER :: i, j, ig, NF, NE, NGW, ldw REAL(DP) :: ANORM REAL(DP) , ALLOCATABLE :: SF(:), SE(:), TEMP(:) COMPLEX(DP), ALLOCATABLE :: CSF(:), CSE(:), CTEMP(:) COMPLEX(DP) :: czero, cone, cmone -! -! ... SUBROUTINE BODY -! - + ! + ! ... SUBROUTINE BODY + ! NF = wfill%nbl( ispin ) NE = wempt%nbl( ispin ) NGW = wfill%ngwl @@ -370,22 +372,24 @@ USE control_flags, ONLY: force_pairing, gamma_only USE wave_functions, ONLY: wave_rand_init -! ... Arguments + ! ... Arguments + COMPLEX(DP), INTENT(INOUT) :: c_occ(:,:,:,:), c_emp(:,:,:,:) TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt REAL(DP) :: ampre -! ... Locals + ! ... Locals + INTEGER :: ig_local INTEGER :: ngw, ngwt INTEGER :: ib, ik, ispin, ispin_wfc LOGICAL :: tortho = .FALSE. COMPLEX(DP), ALLOCATABLE :: pwt( : ) -! -! ... Subroutine body + ! + ! ... Subroutine body -! ... initialize the wave functions in such a way that the values -! ... of the components are independent on the number of processors + ! ... initialize the wave functions in such a way that the values + ! ... of the components are independent on the number of processors ngwt = wfill%ngwt ngw = wfill%ngwl @@ -458,6 +462,8 @@ INTEGER :: i, k, j, iter, ik, nk INTEGER :: nspin, ispin, ispin_wfc INTEGER :: n_occ( wfill%nspin ) + INTEGER :: nupdwn_emp( wfill%nspin ) + INTEGER :: iupdwn_emp( wfill%nspin ) INTEGER :: ig, iprinte, iks, nrl, jl, ngw REAL(DP) :: dek, ekinc, ekinc_old REAL(DP) :: ampre @@ -479,6 +485,10 @@ gamma = wfill%gamma ampre = 0.001d0 + nupdwn_emp = n_emp + iupdwn_emp(1) = 1 + IF( nspin == 2 ) iupdwn_emp(2) = 1+n_emp + ekinc_old = 1.d+10 ekinc = 0.0d0 @@ -517,10 +527,11 @@ CALL nlsm1 ( n_emp, 1, nspnl, eigr, c_emp( 1, 1, ik, ispin ), bece( 1, (ispin-1)*n_emp + 1 ) ) - CALL dforce_all( ispin, c_emp(:,:,1,ispin), wempt, fi(:,1,ispin), eforce(:,:,1,ispin), & - vpot(:,ispin), eigr, bece ) + CALL dforce_all( ispin, c_emp(:,:,1,ispin), fi(:,1,ispin), eforce(:,:,1,ispin), & + vpot(:,ispin), eigr, bece, nupdwn_emp, iupdwn_emp ) ! ... Steepest descent + DO i = 1, n_emp cp_emp(:,i,ik,ispin) = c_emp(:,i,ik,ispin) + dt2bye(:) * eforce(:, i, ik, ispin) END DO @@ -566,7 +577,7 @@ END DO ITERATIONS - CALL empty_eigs( tortho, c_emp, wempt, fi, vpot, eforce, eigr, bece ) + CALL empty_eigs( tortho, c_emp, wempt, fi, vpot, eforce, eigr, bece, nupdwn_emp, iupdwn_emp ) CALL writeempty( c_emp, wempt ) @@ -591,7 +602,7 @@ ! !=----------------------------------------------------------------------------=! - SUBROUTINE empty_eigs( tortho, c_emp, wempt, fi, vpot, eforce, eigr, bece) + SUBROUTINE empty_eigs( tortho, c_emp, wempt, fi, vpot, eforce, eigr, bece, nupdwn_emp, iupdwn_emp) USE wave_types, ONLY : wave_descriptor USE wave_constrains, ONLY : update_lambda @@ -611,6 +622,7 @@ LOGICAL, INTENT(IN) :: TORTHO COMPLEX(DP) :: eigr(:,:) REAL (DP) :: bece(:,:) + INTEGER :: nupdwn_emp(:), iupdwn_emp(:) ! ! ... LOCALS ! @@ -641,8 +653,8 @@ ! ... Calculate | dH / dpsi(j) > ! - CALL dforce_all( ispin, c_emp(:,:,1,ispin), wempt, fi(:,1,ispin), eforce(:,:,1,ispin), & - vpot(:,ispin), eigr, bece ) + CALL dforce_all( ispin, c_emp(:,:,1,ispin), fi(:,1,ispin), eforce(:,:,1,ispin), & + vpot(:,ispin), eigr, bece, nupdwn_emp, iupdwn_emp ) ! ... Calculate Eij = < psi(i) | H | psi(j) > = < psi(i) | dH / dpsi(j) > DO i = 1, n_emp diff --git a/CPV/forces.f90 b/CPV/forces.f90 index bf77bc9d9..664c77450 100644 --- a/CPV/forces.f90 +++ b/CPV/forces.f90 @@ -35,8 +35,8 @@ SUBROUTINE dforce1( co, ce, dco, dce, fio, fie, hg, v, psi_stored ) - USE fft_base, ONLY: dffts - USE gvecw, ONLY: ngw + USE fft_base, ONLY: dffts + USE gvecw, ONLY: ngw USE fft_module, ONLY: fwfft, invfft IMPLICIT NONE @@ -187,7 +187,6 @@ USE ions_base, ONLY: na USE pseudopotential, ONLY: nspnl - USE electrons_base, ONLY: iupdwn USE uspp_param, only: nh USE uspp, only: nhtol, nhtolm, indv, beta, dvan use cvan, only: ish @@ -262,24 +261,20 @@ - SUBROUTINE dforce( ib, iss, c, cdesc, f, df, da, v, eigr, bec ) + SUBROUTINE dforce( ib, iss, c, f, df, da, v, eigr, bec, nupdwn, iupdwn ) ! - USE wave_types, ONLY: wave_descriptor - USE turbo, ONLY: tturbo, nturbo, turbo_states USE reciprocal_vectors, ONLY: ggp, g, gx - USE electrons_base, ONLY: nupdwn, iupdwn ! IMPLICIT NONE ! - INTEGER, INTENT(IN) :: ib, iss ! band and spin index + INTEGER, INTENT(IN) :: ib, iss ! band and spin index COMPLEX(DP), INTENT(IN) :: c(:,:) COMPLEX(DP), INTENT(OUT) :: df(:), da(:) REAL (DP), INTENT(IN) :: v(:), bec(:,:), f(:) COMPLEX(DP), INTENT(IN) :: eigr(:,:) - type (wave_descriptor), INTENT(IN) :: cdesc + INTEGER, INTENT(IN) :: nupdwn(:), iupdwn(:) ! COMPLEX(DP), ALLOCATABLE :: dum( : ) - ! COMPLEX(DP) :: df_( SIZE( df ) ) , da_( SIZE( da ) ) ! DEBUG ! INTEGER :: ig, in ! @@ -295,9 +290,6 @@ ! CALL dforce2_bec( f(ib), f(ib), df , dum , eigr, bec( :, in ), bec( :, in ) ) ! - ! CALL dforce2( f(ib), f(ib), df, dum, fnl(:,:,ib), & - ! fnl(:,:,ib), g(:), gx(:,:), eigr, wsg, wnl ) - ! DEALLOCATE( dum ) ! ELSE @@ -306,14 +298,6 @@ ! CALL dforce2_bec( f(ib), f(ib+1), df, da, eigr, bec( :, in ), bec( :, in+1 ) ) ! - ! CALL dforce2( f(ib), f(ib+1), df, da, fnl(:,:,ib), & - ! fnl(:,:,ib+1), g(:), gx(:,:), eigr, wsg, wnl ) - ! - ! DO ig = 1, SIZE( df ), 50 - ! WRITE(6,*) ig, df(ig), df_(ig) - ! WRITE(6,*) ig, da(ig), da_(ig) - ! END DO - ! END IF ! return @@ -321,34 +305,37 @@ ! ---------------------------------------------- + - - SUBROUTINE dforce_all( ispin, c, cdesc, f, cgrad, vpot, eigr, bec ) - ! - USE wave_types, ONLY: wave_descriptor - USE electrons_base, ONLY: nupdwn + SUBROUTINE dforce_all( iss, c, f, cgrad, vpot, eigr, bec, nupdwn, iupdwn ) ! IMPLICIT NONE - INTEGER, INTENT(IN) :: ispin + INTEGER, INTENT(IN) :: iss COMPLEX(DP), INTENT(INOUT) :: c(:,:) - type (wave_descriptor), INTENT(IN) :: cdesc REAL(DP), INTENT(IN) :: vpot(:), f(:) COMPLEX(DP), INTENT(OUT) :: cgrad(:,:) COMPLEX(DP), INTENT(IN) :: eigr(:,:) REAL(DP), INTENT(IN) :: bec(:,:) + INTEGER, INTENT(IN) :: nupdwn(:), iupdwn(:) INTEGER :: ib ! - IF( nupdwn( ispin ) > 0 ) THEN + IF( nupdwn( iss ) > 0 ) THEN ! ! Process two states at the same time ! - DO ib = 1, nupdwn( ispin ), 2 - CALL dforce( ib, ispin, c, cdesc, f, cgrad(:,ib), cgrad(:,ib+1), & - vpot, eigr, bec ) + DO ib = 1, nupdwn( iss )-1, 2 + CALL dforce( ib, iss, c, f, cgrad(:,ib), cgrad(:,ib+1), & + vpot, eigr, bec, nupdwn, iupdwn ) END DO ! + IF( MOD( nupdwn( iss ), 2 ) /= 0 ) THEN + ib = nupdwn( iss ) + CALL dforce( ib, iss, c, f, cgrad(:,ib), cgrad(:,ib), & + vpot, eigr, bec, nupdwn, iupdwn ) + END IF + ! END IF ! RETURN diff --git a/CPV/ksstates.f90 b/CPV/ksstates.f90 index e52ef90a7..f055c7f4f 100644 --- a/CPV/ksstates.f90 +++ b/CPV/ksstates.f90 @@ -163,6 +163,7 @@ CONTAINS USE brillouin, ONLY: kpoints, kp USE pseudo_projector, ONLY: projector USE control_flags, ONLY: timing, force_pairing + USE electrons_base, ONLY: nupdwn, iupdwn IMPLICIT NONE @@ -176,6 +177,8 @@ CONTAINS ! ... declare other variables INTEGER :: i, ik, ib, nk, ig, ngw, nb_g, nb_l, ispin, nspin, iks INTEGER :: ispin_wfc + INTEGER :: nupdwn_emp( wfill%nspin ) + INTEGER :: iupdwn_emp( wfill%nspin ) LOGICAL :: tortho = .TRUE. CHARACTER(LEN=4) :: nom CHARACTER(LEN=256) :: file_name @@ -211,8 +214,8 @@ CONTAINS ALLOCATE( eforce( ngw, nb_l ) ) - CALL dforce_all( ispin, cf(:,:,1,ispin_wfc), wfill, occ(:,1,ispin), eforce, & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, cf(:,:,1,ispin_wfc), occ(:,1,ispin), eforce, & + vpot(:,ispin), eigr, bec, nupdwn, iupdwn ) CALL kohn_sham( ispin, cf(:,:,1,ispin_wfc), wfill, eforce ) @@ -225,6 +228,10 @@ CONTAINS ngw = wempt%ngwl nb_l = wempt%nbl( ispin ) + nupdwn_emp = nb_l + iupdwn_emp(1) = 1 + IF( nspin == 2 ) iupdwn_emp(2) = 1+nb_l + IF( nb_l > 0 ) THEN ALLOCATE( fi( nb_l, nk ) ) @@ -234,8 +241,8 @@ CONTAINS ALLOCATE( eforce( ngw, nb_l ) ) - CALL dforce_all( ispin, ce(:,:,1,ispin), wempt, fi(:,1), eforce, & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, ce(:,:,1,ispin), fi(:,1), eforce, & + vpot(:,ispin), eigr, bec, nupdwn_emp, iupdwn_emp ) CALL kohn_sham( ispin, ce(:,:,1,ispin), wempt, eforce ) @@ -357,7 +364,7 @@ CONTAINS USE brillouin, ONLY: kpoints, kp USE pseudo_projector, ONLY: projector USE control_flags, ONLY: timing - USE electrons_module, ONLY: nupdwn, nspin + USE electrons_base, ONLY: iupdwn, nupdwn, nspin IMPLICIT NONE @@ -378,6 +385,8 @@ CONTAINS COMPLEX(DP), ALLOCATABLE :: eforce(:,:,:) REAL(DP), ALLOCATABLE :: fi(:,:) + INTEGER :: nupdwn_emp( wfill%nspin ) + INTEGER :: iupdwn_emp( wfill%nspin ) CHARACTER (LEN=6), EXTERNAL :: int_to_char REAL(DP), EXTERNAL :: cclock @@ -406,10 +415,8 @@ CONTAINS ALLOCATE( eforce( ngw, nb, 2 ) ) - CALL dforce_all( 1, cf(:,:,1,1), wfill, occ(:,1,1), eforce(:,:,1), & - vpot(:,1), eigr, bec ) - CALL dforce_all( 2, cf(:,:,1,1), wfill, occ(:,1,2), eforce(:,:,2), & - vpot(:,2), eigr, bec ) + CALL dforce_all( 1, cf(:,:,1,1), occ(:,1,1), eforce(:,:,1), vpot(:,1), eigr, bec, nupdwn, iupdwn ) + CALL dforce_all( 2, cf(:,:,1,1), occ(:,1,2), eforce(:,:,2), vpot(:,2), eigr, bec, nupdwn, iupdwn ) DO i = 1, nupdwn(2) eforce(:,i,1) = occ(i,1,1) * eforce(:,i,1) + occ(i,1,2) * eforce(:,i,2) @@ -429,6 +436,10 @@ CONTAINS ngw = wempt%ngwl nb_l = wempt%nbl( 1 ) + nupdwn_emp = nb_l + iupdwn_emp(1) = 1 + IF( nspin == 2 ) iupdwn_emp(2) = 1+nb_l + IF( nb_l > 0 ) THEN ALLOCATE( fi( nb_l, nk ) ) @@ -438,13 +449,11 @@ CONTAINS ALLOCATE( eforce( ngw, nb_l, 1 )) - CALL dforce_all( 1, ce(:,:,1,1), wempt, fi(:,1), eforce(:,:,1), vpot(:,1), & - eigr, bec ) + CALL dforce_all( 1, ce(:,:,1,1), fi(:,1), eforce(:,:,1), vpot(:,1), eigr, bec, nupdwn_emp, iupdwn_emp ) CALL kohn_sham( 1, ce(:,:,1,1), wempt, eforce(:,:,1) ) - CALL dforce_all( 2, ce(:,:,1,2), wempt, fi(:,1), eforce(:,:,1), vpot(:,2), & - eigr, bec ) + CALL dforce_all( 2, ce(:,:,1,2), fi(:,1), eforce(:,:,1), vpot(:,2), eigr, bec, nupdwn_emp, iupdwn_emp ) CALL kohn_sham( 2, ce(:,:,1,2), wempt, eforce(:,:,1) ) diff --git a/CPV/optical.f90 b/CPV/optical.f90 index ad9245f96..71d606a79 100644 --- a/CPV/optical.f90 +++ b/CPV/optical.f90 @@ -68,6 +68,7 @@ USE forces, ONLY: dforce_all USE brillouin, ONLY: kpoints, kp USE electrons_module, ONLY: ei, ei_emp + USE electrons_base, ONLY: nupdwn, iupdwn USE kohn_sham_states, ONLY: kohn_sham USE constants, ONLY: au, pi, k_boltzman_au, au_to_ohmcmm1 USE cell_base, ONLY: tpiba2 @@ -107,6 +108,10 @@ REAL(DP) :: FACT, ef, beta LOGICAL :: gamma_symmetry, gzero + INTEGER :: nupdwn_emp( wfill%nspin ) + INTEGER :: iupdwn_emp( wfill%nspin ) + + ! ... SUBROUTINE BODY CALL errore( ' opticalp ', ' not working ', 1 ) @@ -138,13 +143,17 @@ ! CALL nlsm1 ( nb_l, 1, nspnl, eigr, cf(1,1,1,ispin), bec ) - CALL dforce_all( ispin, cf(:,:,1,ispin), wfill, occ(:,1,ispin), eforce(:,:,1), & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, cf(:,:,1,ispin), occ(:,1,ispin), eforce(:,:,1), & + vpot(:,ispin), eigr, bec, nupdwn, iupdwn ) CALL kohn_sham( ispin, cf(:,:,1,ispin), wfill, eforce(:,:,1) ) ngw = wempt%ngwl nb_l = wempt%nbl( ispin ) + nupdwn_emp = nb_l + iupdwn_emp(1) = 1 + IF( nspin == 2 ) iupdwn_emp(2) = 1+nb_l + ALLOCATE( ff( nb_l, nk ) ) DO ik = 1, nk @@ -155,8 +164,8 @@ ! CALL nlsm1 ( nb_l, 1, nspnl, eigr, ce(1,1,1,ispin), bece ) ! - CALL dforce_all( ispin, ce(:,:,1,ispin), wempt, ff( :, 1), eforce(:,:,1), & - vpot(:,ispin), eigr, bece ) + CALL dforce_all( ispin, ce(:,:,1,ispin), ff( :, 1), eforce(:,:,1), & + vpot(:,ispin), eigr, bece, nupdwn_emp, iupdwn_emp ) ! CALL kohn_sham( ispin, ce(:,:,1,ispin), wempt, eforce(:,:,1) ) diff --git a/CPV/runcg.f90 b/CPV/runcg.f90 index 68675c347..24ef62ee0 100644 --- a/CPV/runcg.f90 +++ b/CPV/runcg.f90 @@ -63,6 +63,7 @@ ! ... declare modules USE energies, ONLY: dft_energy_type, print_energies USE electrons_module, ONLY: pmss, eigs, nb_l + USE electrons_base, ONLY: nupdwn, iupdwn USE cp_electronic_mass, ONLY: emass USE wave_functions, ONLY: cp_kinetic_energy, proj, fixwave USE wave_base, ONLY: dotp, hpsi @@ -183,8 +184,8 @@ ! ... Calculate wave functions gradient (temporarely stored in cp) ! ... |d H / dPsi_j > = H |Psi_j> - Sum{i} |Psi_i> - CALL dforce_all( ispin, c0(:,:,1,ispin), cdesc, occ(:,1,ispin), cp(:,:,1,ispin), & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, c0(:,:,1,ispin), occ(:,1,ispin), cp(:,:,1,ispin), & + vpot(:,ispin), eigr, bec, nupdwn, iupdwn ) ! ... Project the gradient IF( gamma_symmetry ) THEN @@ -317,8 +318,8 @@ IF( tprint ) THEN DO ispin = 1, nspin - CALL dforce_all( ispin, c0(:,:,1,ispin), cdesc, occ(:,1,ispin), hacca(:,:,1,ispin), & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, c0(:,:,1,ispin), occ(:,1,ispin), hacca(:,:,1,ispin), & + vpot(:,ispin), eigr, bec, nupdwn, iupdwn ) nb_g( ispin ) = cdesc%nbt( ispin ) diff --git a/CPV/runcp.f90 b/CPV/runcp.f90 index 1ca69770c..a77150822 100644 --- a/CPV/runcp.f90 +++ b/CPV/runcp.f90 @@ -47,7 +47,6 @@ USE wave_base, ONLY: hpsi USE cell_module, ONLY: boxdimensions USE time_step, ONLY: delt - USE forces, ONLY: dforce USE orthogonalize, ONLY: ortho USE wave_types, ONLY: wave_descriptor USE pseudo_projector, ONLY: projector @@ -167,6 +166,7 @@ ! ... declare modules USE kinds USE electrons_module, ONLY: pmss + USE electrons_base, ONLY: nupdwn, iupdwn USE cp_electronic_mass, ONLY: emass USE wave_base, ONLY: wave_steepest, wave_verlet USE time_step, ONLY: delt @@ -252,7 +252,7 @@ DO i = 1, nb, 2 !WRITE( 6, * ) 'DEBUG = ', fi(i,1,is), fi(i+1,1,is) - CALL dforce( i, is, c0(:,:,1,is), cdesc, fi(:,1,is), c2, c3, vpot(:,is), eigr, bec ) + CALL dforce( i, is, c0(:,:,1,is), fi(:,1,is), c2, c3, vpot(:,is), eigr, bec, nupdwn, iupdwn ) IF( tlam ) THEN CALL update_lambda( i, gam( :, :,is), c0(:,:,1,is), cdesc, c2 ) @@ -283,7 +283,7 @@ nb = nx - CALL dforce( nx, is, c0(:,:,1,is), cdesc, fi(:,1,is), c2, c3, vpot(:,is), eigr, bec ) + CALL dforce( nx, is, c0(:,:,1,is), fi(:,1,is), c2, c3, vpot(:,is), eigr, bec, nupdwn, iupdwn ) IF( tlam ) THEN CALL update_lambda( nb, gam( :, :,is), c0(:,:,1,is), cdesc, c2 ) @@ -332,7 +332,8 @@ USE kinds USE mp_global, ONLY: intra_image_comm USE mp, ONLY: mp_sum - USE electrons_module, ONLY: pmss, eigs, nb_l, nupdwn, nspin + USE electrons_module, ONLY: pmss, eigs, nb_l + USE electrons_base, ONLY: iupdwn, nupdwn, nspin USE cp_electronic_mass, ONLY: emass USE wave_functions, ONLY : cp_kinetic_energy USE wave_base, ONLY: wave_steepest, wave_verlet @@ -454,8 +455,8 @@ DO i = 1, nb, 2 ! - CALL dforce( i, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,1), eigr, bec ) - CALL dforce( i, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c4, c5, vpot(:,2), eigr, bec ) + CALL dforce( i, 2, c0(:,:,1,1), fi(:,1,1), c2, c3, vpot(:,1), eigr, bec, nupdwn, iupdwn ) + CALL dforce( i, 2, c0(:,:,1,1), fi(:,1,1), c4, c5, vpot(:,2), eigr, bec, nupdwn, iupdwn ) ! c2 = occup(i , ik)* (c2 + c4) c3 = occup(i+1, ik)* (c3 + c5) @@ -487,8 +488,8 @@ ! nb = n_unp - 1 ! - CALL dforce( nb, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,1), eigr, bec ) - CALL dforce( nb, 2, c0(:,:,1,1), cdesc, fi(:,1,2), c4, c5, vpot(:,2), eigr, bec ) + CALL dforce( nb, 2, c0(:,:,1,1), fi(:,1,1), c2, c3, vpot(:,1), eigr, bec, nupdwn, iupdwn ) + CALL dforce( nb, 2, c0(:,:,1,1), fi(:,1,2), c4, c5, vpot(:,2), eigr, bec, nupdwn, iupdwn ) c2 = occup(nb , ik)* (c2 + c4) @@ -506,7 +507,7 @@ END IF ! - CALL dforce( n_unp, 1, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,1), eigr, bec ) + CALL dforce( n_unp, 1, c0(:,:,1,1), fi(:,1,1), c2, c3, vpot(:,1), eigr, bec, nupdwn, iupdwn ) intermed = -2.d0 * sum( c2 * conjg( c0(:, n_unp, ik, 1 ) ) ) IF ( gstart == 2 ) THEN diff --git a/CPV/rundiis.f90 b/CPV/rundiis.f90 index 6ba86f1c8..63fb905ed 100644 --- a/CPV/rundiis.f90 +++ b/CPV/rundiis.f90 @@ -79,6 +79,7 @@ CONTAINS USE io_global, ONLY: stdout USE energies, ONLY: dft_energy_type USE electrons_module, ONLY: pmss + USE electrons_base, ONLY: nupdwn, iupdwn USE time_step, ONLY: delt USE wave_functions, ONLY: proj, crot USE phase_factors_module, ONLY: strucf, phfacs @@ -263,7 +264,7 @@ CONTAINS edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, bec, becdr, eigr) - CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec ) + CALL dforce_all( 1, c0(:,:,1,1), fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec, nupdwn, iupdwn ) CALL proj( 1, cgrad(:,:,1,1), cdesc, c0(:,:,1,1), cdesc, lambda ) CALL crot( 1, c0(:,:,1,1), cdesc, lambda, eig(:,1,1) ) @@ -272,7 +273,7 @@ CONTAINS call entropy_s(fi(1,1,1),temp_elec,cdesc%nbl(1),edft%ent) edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, bec, becdr, eigr) - CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec ) + CALL dforce_all( 1, c0(:,:,1,1), fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec, nupdwn, iupdwn ) DO ib = 1, cdesc%nbl( 1 ) cgrad(:,ib,1,1) = cgrad(:,ib,1,1) + eig(ib,1,1)*c0(:,ib,1,1) @@ -283,7 +284,7 @@ CONTAINS ! ... DIIS on c0 at FIXED potential edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, bec, becdr, eigr) - CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec ) + CALL dforce_all( 1, c0(:,:,1,1), fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec, nupdwn, iupdwn ) CALL proj( 1, cgrad(:,:,1,1), cdesc, c0(:,:,1,1), cdesc, lambda) @@ -393,6 +394,7 @@ CONTAINS USE runcp_module, ONLY: runcp USE energies, ONLY: dft_energy_type USE electrons_module, ONLY: ei, pmss + USE electrons_base, ONLY: nupdwn, iupdwn USE time_step, ONLY: delt USE wave_functions, ONLY: proj, update_wave_functions USE diis @@ -528,8 +530,8 @@ CONTAINS ! ... of 1, ( row 1 to PE 1, row 2 to PE 2, .. row nproc_image+1 to PE 1 and ! ... so on). - CALL dforce_all( ispin, c0(:,:,1,ispin), cdesc, fi(:,1,ispin), cgrad(:,:,1,ispin), & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, c0(:,:,1,ispin), fi(:,1,ispin), cgrad(:,:,1,ispin), & + vpot(:,ispin), eigr, bec, nupdwn, iupdwn ) CALL proj( ispin, cgrad(:,:,1,ispin), cdesc, c0(:,:,1,ispin), cdesc, lambda) @@ -604,6 +606,7 @@ CONTAINS USE constants, ONLY: au USE cell_base, ONLY: tpiba2 USE electrons_module, ONLY: eigs, ei, pmss, emass, nb_l, ib_owner, ib_local + USE electrons_base, ONLY: iupdwn, nupdwn USE forces, ONLY: dforce_all USE orthogonalize USE pseudopotential, ONLY: nspnl @@ -653,8 +656,8 @@ CONTAINS CALL nlsm1( n, 1, nspnl, eigr, c(1,1,1,ispin), bec ) ! ... Calculate | dH / dpsi(j) > - CALL dforce_all( ispin, c(:,:,1,ispin), cdesc, fi(:,1,ispin), eforce(:,:,1,ispin), & - vpot(:,ispin), eigr, bec ) + CALL dforce_all( ispin, c(:,:,1,ispin), fi(:,1,ispin), eforce(:,:,1,ispin), & + vpot(:,ispin), eigr, bec, nupdwn, iupdwn ) ! ... Calculate Eij = < psi(i) | H | psi(j) > = < psi(i) | dH / dpsi(j) > DO i = 1, n