mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-29 06:35:28 -04:00
Clean up. Fixed few bugs with parallel. Chemical shifts are now computed by adding up g-coeffs (seems more accurate). Interpolation based method to compute the shifts needs to be carefuly parallelized (problem with the distributed grid).
svn-origin-rev: 6085
This commit is contained in:
Valery Weber
parent
53dc0da0ce
commit
52a0a5002a
4 changed files with 464 additions and 399 deletions
|
|
@ -110,6 +110,7 @@ MODULE qs_linres_nmr_current
|
|||
set_vecp_rev
|
||||
USE qs_linres_nmr_shift, ONLY: chi_soft_analytic,&
|
||||
interpolate_shift_pwgrid,&
|
||||
gsum_shift_pwgrid,&
|
||||
mult_G_ov_G2_grid
|
||||
USE qs_linres_types, ONLY: nmr_env_type
|
||||
USE qs_mo_types, ONLY: get_mo_set,&
|
||||
|
|
@ -194,22 +195,19 @@ CONTAINS
|
|||
!!*** **********************************************************************
|
||||
|
||||
SUBROUTINE nmr_response_current(nmr_env,qs_env,psi1,p_psi1,error)
|
||||
|
||||
!
|
||||
TYPE(nmr_env_type) :: nmr_env
|
||||
TYPE(qs_environment_type), POINTER :: qs_env
|
||||
TYPE(cp_fm_p_type), DIMENSION(:), &
|
||||
POINTER :: psi1, p_psi1
|
||||
POINTER :: psi1, p_psi1
|
||||
TYPE(cp_error_type), INTENT(INOUT) :: error
|
||||
|
||||
!
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'nmr_response_current', &
|
||||
routineP = moduleN//':'//routineN
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
CHARACTER(LEN=80) :: ext, filename
|
||||
CHARACTER, DIMENSION(3) :: labels
|
||||
|
||||
DATA labels /'x','y','z'/
|
||||
INTEGER :: handle, homo, i_B, iao, iatom, idir, idir2, idir3, ii_B, &
|
||||
iii_B, ispin, istat, istate, nao, natom, nmo, nspins, unit_nr, output_unit
|
||||
iii_B, ispin, istat, istate, nao, natom, nmo, nspins, unit_nr, output_unit
|
||||
LOGICAL :: failure, gapw, ionode, uni_occ
|
||||
REAL(dp) :: alpha, dk(3), jrho_tot, &
|
||||
maxocc, my_chi,&
|
||||
|
|
@ -248,10 +246,10 @@ CONTAINS
|
|||
TYPE(realspace_grid_type), POINTER :: rs
|
||||
TYPE(section_vals_type), POINTER :: nmr_section
|
||||
REAL(dp), EXTERNAL :: DDOT
|
||||
! ---------------------------------------------------------------------------
|
||||
|
||||
!
|
||||
failure = .FALSE.
|
||||
gapw=.FALSE.
|
||||
CALL timeset("nmr_response_current","I"," ",handle)
|
||||
!
|
||||
NULLIFY(rho_rspace,rho_gspace,logger,nmr_section)
|
||||
NULLIFY(density_matrix,density_matrix_ii,density_matrix_iii)
|
||||
NULLIFY(cell,dft_control,mos,rho_struct,particle_set,pw_env)
|
||||
|
|
@ -266,374 +264,364 @@ CONTAINS
|
|||
nspins = SIZE(psi1,1)
|
||||
ALLOCATE(density_matrix(nspins),density_matrix_ii(nspins),density_matrix_iii(nspins))
|
||||
DO ispin = 1,nspins
|
||||
NULLIFY(density_matrix(ispin)%matrix)
|
||||
NULLIFY(density_matrix_ii(ispin)%matrix)
|
||||
NULLIFY(density_matrix_iii(ispin)%matrix)
|
||||
END DO
|
||||
|
||||
|
||||
NULLIFY(density_matrix(ispin)%matrix)
|
||||
NULLIFY(density_matrix_ii(ispin)%matrix)
|
||||
NULLIFY(density_matrix_iii(ispin)%matrix)
|
||||
ENDDO
|
||||
!
|
||||
CALL get_qs_env(qs_env=qs_env,rho=rho_struct,&
|
||||
cell=cell, dft_control=dft_control,mos=mos,&
|
||||
para_env=para_env,particle_set=particle_set,error=error)
|
||||
gapw = dft_control%qs_control%gapw
|
||||
& cell=cell,dft_control=dft_control,mos=mos,&
|
||||
& para_env=para_env,particle_set=particle_set,error=error)
|
||||
gapw = dft_control%qs_control%gapw
|
||||
natom = SIZE(particle_set,1)
|
||||
CPPrecondition(ASSOCIATED(rho_struct),cp_failure_level,routineP,error,failure)
|
||||
CPPrecondition(rho_struct%ref_count>0,cp_failure_level,routineP,error,failure)
|
||||
nmr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES%NMR",&
|
||||
error=error)
|
||||
|
||||
& error=error)
|
||||
!
|
||||
IF (.NOT. failure) THEN
|
||||
CALL timeset("nmr_response_current","I"," ",handle)
|
||||
|
||||
CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error)
|
||||
CALL pw_env_get(pw_env, auxbas_rs_pool=auxbas_rs_pool,&
|
||||
auxbas_pw_pool=auxbas_pw_pool,pw_pools=pw_pools,error=error)
|
||||
IF (BTEST(cp_print_key_should_output(logger%iter_info,nmr_section,&
|
||||
"PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN
|
||||
|
||||
CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error)
|
||||
CALL pw_env_get(pw_env, auxbas_rs_pool=auxbas_rs_pool,&
|
||||
& auxbas_pw_pool=auxbas_pw_pool,pw_pools=pw_pools,error=error)
|
||||
IF (BTEST(cp_print_key_should_output(logger%iter_info,nmr_section,&
|
||||
& "PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN
|
||||
CALL rs_pool_create_rs(auxbas_rs_pool,rs,error=error)
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,wf_r,use_data=REALDATA3D,&
|
||||
& in_space=REALSPACE,error=error)
|
||||
ENDIF
|
||||
!
|
||||
! Allocate grids for the calculation of jrho and the shift
|
||||
ALLOCATE(shift_pw_gspace(3),STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
DO idir = 1,3
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,shift_pw_gspace(idir),&
|
||||
& use_data=COMPLEXDATA1D,&
|
||||
& in_space=RECIPROCALSPACE,error=error)
|
||||
CALL coeff_zero(shift_pw_gspace(idir))
|
||||
ENDDO
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,shift_pw_rspace,&
|
||||
& use_data=REALDATA3D,in_space=REALSPACE,error=error)
|
||||
CALL coeff_zero(shift_pw_rspace)
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,pw_gspace_work,&
|
||||
& use_data=COMPLEXDATA1D,&
|
||||
& in_space=RECIPROCALSPACE,error=error)
|
||||
CALL coeff_zero(pw_gspace_work)
|
||||
!
|
||||
! Initialize shift and chi
|
||||
nmr_env%chemical_shift = 0.0_dp
|
||||
nmr_env%chemical_shift_loc = 0.0_dp
|
||||
IF(nmr_env%do_nics) THEN
|
||||
nmr_env%chemical_shift_nics = 0.0_dp
|
||||
nmr_env%chemical_shift_loc_nics = 0.0_dp
|
||||
ENDIF
|
||||
nmr_env%chi_tensor = 0.0_dp
|
||||
nmr_env%chi_tensor_loc = 0.0_dp
|
||||
jrho_tot_G = 0.0_dp
|
||||
jrho_tot_R = 0.0_dp
|
||||
!
|
||||
! Loop on the field direction
|
||||
DO i_B = 1,3
|
||||
DO ispin = 1,dft_control%nspins
|
||||
CALL cp_fm_set_all(psi1(ispin)%matrix,0.0_dp,error=error)
|
||||
CALL cp_fm_set_all(p_psi1(ispin)%matrix,0.0_dp,error=error)
|
||||
ENDDO
|
||||
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,wf_r,&
|
||||
use_data = REALDATA3D,&
|
||||
in_space = REALSPACE, error=error)
|
||||
END IF
|
||||
CALL set_vecp(i_B,ii_B,iii_B)
|
||||
DO ispin = 1,nspins
|
||||
CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff,&
|
||||
& occupation_numbers=occupation,homo=homo,nao=nao,nmo=nmo,&
|
||||
& uniform_occupation=uni_occ,maxocc=maxocc)
|
||||
NULLIFY(fm_work1)
|
||||
!
|
||||
! create a new matrix
|
||||
NULLIFY(tmp_fm_struct)
|
||||
CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,&
|
||||
& ncol_global=nmo,para_env=para_env,&
|
||||
& context=mo_coeff%matrix_struct%context,&
|
||||
& error=error)
|
||||
IF(ASSOCIATED(fm_work1))THEN
|
||||
CALL cp_fm_release(fm_work1,error=error)
|
||||
ENDIF
|
||||
CALL cp_fm_create(fm_work1,tmp_fm_struct,error=error)
|
||||
CALL cp_fm_set_all(fm_work1,0.0_dp,error=error)
|
||||
IF ( .NOT. uni_occ ) THEN
|
||||
CALL cp_fm_create(psi0_fk,tmp_fm_struct,error=error)
|
||||
CALL cp_fm_set_all(psi0_fk,0.0_dp,error=error)
|
||||
ELSE
|
||||
NULLIFY(psi0_fk)
|
||||
END IF
|
||||
CALL cp_fm_struct_release(tmp_fm_struct,error=error)
|
||||
!
|
||||
! Allocate buffer vectors
|
||||
DO idir = 1,3
|
||||
ALLOCATE(vecbuf(idir)%array(1,nao),STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
END DO
|
||||
ALLOCATE(vecbuf_psi1(1,nao),STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
!
|
||||
! First compute the contribution to the CHI tensor coming from the soft density
|
||||
CALL chi_soft_analytic(nmr_env,qs_env,ispin,i_B,error=error)
|
||||
!
|
||||
! Construct the 3 density matrices for the field in direction i_B
|
||||
!
|
||||
! First the full matrix psi_a_iB
|
||||
! psi1_rxp - psi1_D - (R_\nu-dk)_ii psi1_piiiB + (R_\nu-dk)_iii psi1_piiB
|
||||
psi_a_iB => psi1(ispin)%matrix
|
||||
CALL cp_fm_set_all(psi_a_iB,0.0_dp,error=error)
|
||||
!
|
||||
! contributions from the response psi1_p_ii and psi1_p_iii
|
||||
DO istate = 1,nmr_env%nstates(ispin)
|
||||
CALL cp_fm_get_submatrix(nmr_env%psi1_p(ispin,ii_B)%matrix,&
|
||||
& vecbuf(ii_B)%array,1,istate,nao,1,&
|
||||
& transpose=.TRUE.,error=error)
|
||||
CALL cp_fm_get_submatrix(nmr_env%psi1_p(ispin,iii_B)%matrix,&
|
||||
& vecbuf(iii_B)%array,1,istate,nao,1,&
|
||||
& transpose=.TRUE.,error=error)
|
||||
dk(1:3) = nmr_env%centers_set(ispin)%array(1:3,istate)
|
||||
DO iao = 1,nao
|
||||
Rmu(1:3) = nmr_env%basisfun_center(1:3,iao)
|
||||
Rmu_dk = pbc(dk,Rmu,cell)
|
||||
vecbuf_psi1(1,iao) = Rmu_dk(ii_B)*vecbuf(iii_B)%array(1,iao)-&
|
||||
& Rmu_dk(iii_B)*vecbuf(ii_B)%array(1,iao)
|
||||
ENDDO ! iao
|
||||
!
|
||||
! Copy the vector in the full matrix psi1
|
||||
CALL cp_fm_set_submatrix(psi1(ispin)%matrix,vecbuf_psi1,&
|
||||
& 1,istate,nao,1,transpose=.TRUE.,&
|
||||
& error=error)
|
||||
ENDDO ! istate
|
||||
!
|
||||
! contribution from the response psi1_rxp
|
||||
CALL cp_fm_scale_and_add(-1.0_dp,psi_a_iB,1.0_dp,&
|
||||
& nmr_env%psi1_rxp(ispin,i_B)%matrix,&
|
||||
& error=error)
|
||||
IF(nmr_env%full_nmr) THEN
|
||||
! contribution from the response psi1_D
|
||||
CALL cp_fm_scale_and_add(1.0_dp,psi_a_iB,-1.0_dp,&
|
||||
& nmr_env%psi1_D(ispin,i_B)%matrix,&
|
||||
& error=error)
|
||||
ENDIF
|
||||
!
|
||||
! Multiply by the occupation number for the density matrix
|
||||
IF ( .NOT. uni_occ ) THEN
|
||||
alpha = 1.0_dp
|
||||
CALL cp_fm_to_fm(mo_coeff,psi0_fk,error=error)
|
||||
CALL cp_fm_column_scale(psi0_fk,occupation(1:homo))
|
||||
ELSE
|
||||
alpha = maxocc
|
||||
psi0_fk => mo_coeff
|
||||
ENDIF
|
||||
!
|
||||
! Build the first density matrix
|
||||
density_matrix(ispin)%matrix => nmr_env%jp1_ao
|
||||
CALL set_matrix(density_matrix(ispin)%matrix ,0.0_dp)
|
||||
CALL cp_sm_plus_fm_fm_t(sparse_matrix=density_matrix(ispin)%matrix ,&
|
||||
& matrix_v=psi0_fk,matrix_g=psi_a_iB,&
|
||||
& ncol=homo,alpha=alpha,error=error)
|
||||
!
|
||||
! Build the second density matrix
|
||||
density_matrix_iii(ispin)%matrix => nmr_env%jp2_ao(1)%matrix
|
||||
CALL set_matrix(density_matrix_iii(ispin)%matrix,0.0_dp)
|
||||
CALL cp_sm_plus_fm_fm_t(sparse_matrix=density_matrix_iii(ispin)%matrix ,&
|
||||
& matrix_v=psi0_fk,matrix_g=nmr_env%psi1_p(ispin,iii_B)%matrix,&
|
||||
& ncol=homo,alpha=alpha,error=error)
|
||||
!
|
||||
! Build the third density matrix
|
||||
density_matrix_ii(ispin)%matrix => nmr_env%jp2_ao(2)%matrix
|
||||
CALL set_matrix(density_matrix_ii(ispin)%matrix ,0.0_dp)
|
||||
CALL cp_sm_plus_fm_fm_t(sparse_matrix=density_matrix_ii(ispin)%matrix ,&
|
||||
& matrix_v=psi0_fk,matrix_g=nmr_env%psi1_p(ispin,ii_B)%matrix,&
|
||||
& ncol=homo,alpha=alpha,error=error)
|
||||
DO idir3=1,3
|
||||
! set to zero for the calculation of the shift
|
||||
CALL coeff_zero(shift_pw_gspace(idir3))
|
||||
ENDDO
|
||||
DO idir = 1,3
|
||||
!
|
||||
! Calculate the current density on the pw grid (only soft if GAPW)
|
||||
! idir is the cartesian component of the response current density
|
||||
! generated by the magnetic field pointing in cartesian direction i_B
|
||||
! Use the qs_rho_type already used for rho during the scf
|
||||
|
||||
IF(nmr_env%store_current) THEN
|
||||
rho_rspace => nmr_env%jrho1_set(idir,i_B)%rho%rho_r(ispin)
|
||||
rho_gspace => nmr_env%jrho1_set(idir,i_B)%rho%rho_g(ispin)
|
||||
ELSE
|
||||
rho_rspace => rho_struct%rho_r(ispin)
|
||||
rho_gspace => rho_struct%rho_g(ispin)
|
||||
ENDIF
|
||||
|
||||
!-----------------------------------------------------------------------!
|
||||
! Allocate grids for the calculation of jrho and the shift
|
||||
ALLOCATE(shift_pw_gspace(3),STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
DO idir = 1,3
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,shift_pw_gspace(idir), &
|
||||
use_data = COMPLEXDATA1D,&
|
||||
in_space = RECIPROCALSPACE, error=error)
|
||||
CALL coeff_zero(shift_pw_gspace(idir))
|
||||
END DO
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,shift_pw_rspace,&
|
||||
use_data=REALDATA3D, in_space=REALSPACE, error=error)
|
||||
CALL coeff_zero(shift_pw_rspace)
|
||||
CALL pw_pool_init_coeff(auxbas_pw_pool,pw_gspace_work,&
|
||||
use_data = COMPLEXDATA1D,&
|
||||
in_space = RECIPROCALSPACE,error=error)
|
||||
CALL coeff_zero(pw_gspace_work)
|
||||
!
|
||||
!-----------------------------------------------------------------------!
|
||||
CALL coeff_zero(rho_rspace)
|
||||
CALL coeff_zero(rho_gspace)
|
||||
CALL calculate_jrho_resp(density_matrix(ispin)%matrix, &
|
||||
& density_matrix_ii(ispin)%matrix,&
|
||||
& density_matrix_iii(ispin)%matrix, &
|
||||
& i_B,idir,rho_rspace,rho_gspace,qs_env,&
|
||||
& gapw,error=error)
|
||||
|
||||
scale_fac = cell%deth / twopi
|
||||
CALL coeff_scale(rho_rspace,scale_fac)
|
||||
CALL coeff_scale(rho_gspace,scale_fac)
|
||||
|
||||
!-----------------------------------------------------------------------!
|
||||
! Initialize shift and chi
|
||||
nmr_env%chemical_shift = 0.0_dp
|
||||
nmr_env%chemical_shift_loc = 0.0_dp
|
||||
IF(nmr_env%do_nics) THEN
|
||||
nmr_env%chemical_shift_nics = 0.0_dp
|
||||
nmr_env%chemical_shift_loc_nics = 0.0_dp
|
||||
END IF
|
||||
nmr_env%chi_tensor = 0.0_dp
|
||||
nmr_env%chi_tensor_loc = 0.0_dp
|
||||
jrho_tot_G = 0.0_dp
|
||||
jrho_tot_R = 0.0_dp
|
||||
!-----------------------------------------------------------------------!
|
||||
jrho_tot_G(idir,i_B) = calculate_total_rho(rho_gspace)
|
||||
jrho_tot_R(idir,i_B) = calculate_total_rho(rho_rspace)
|
||||
|
||||
! Loop on the field direction
|
||||
DO i_B = 1,3
|
||||
DO ispin = 1,dft_control%nspins
|
||||
CALL cp_fm_set_all(psi1(ispin)%matrix,0.0_dp,error=error)
|
||||
CALL cp_fm_set_all(p_psi1(ispin)%matrix,0.0_dp,error=error)
|
||||
END DO
|
||||
IF(output_unit>0) THEN
|
||||
WRITE(output_unit,'(T2,2(A,E24.16))') 'Integrated j_'&
|
||||
&//ACHAR(idir+119)//ACHAR(i_B+119)//'(r): G-space=',&
|
||||
&jrho_tot_G(idir,i_B),' R-space=',jrho_tot_R(idir,i_B)
|
||||
ENDIF
|
||||
!
|
||||
! Field gradient
|
||||
! loop over the Gvec components: x,y,z
|
||||
DO idir2 = 1,3
|
||||
IF(idir /= idir2) THEN
|
||||
! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i))
|
||||
! We dont need to set my_chi here, it is done while printing
|
||||
! the chemical shifts.
|
||||
my_chi = 0.0_dp !nmr_env%chi_SI2shiftppm/nmr_env%shift_factor*&
|
||||
! nmr_env%chi_tensor(idir3,i_B)
|
||||
CALL mult_G_ov_G2_grid(cell,auxbas_pw_pool,rho_gspace,&
|
||||
& pw_gspace_work,idir2,my_chi,error=error)
|
||||
!
|
||||
! scale and add to the correct component of the shift column
|
||||
CALL set_vecp_rev(idir,idir2,idir3)
|
||||
scale_fac=fac_vecp(idir3,idir2,idir)
|
||||
CALL coeff_scale(pw_gspace_work,scale_fac)
|
||||
CALL coeff_sumup(pw_gspace_work,shift_pw_gspace(idir3))
|
||||
ENDIF
|
||||
ENDDO
|
||||
ENDDO ! idir
|
||||
!
|
||||
! Deallocate buffer vectors
|
||||
DO idir = 1,3
|
||||
DEALLOCATE(vecbuf(idir)%array,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
ENDDO
|
||||
DEALLOCATE(vecbuf_psi1,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
CALL set_vecp(i_B,ii_B,iii_B)
|
||||
DO ispin = 1,nspins
|
||||
CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff,&
|
||||
occupation_numbers=occupation, homo=homo, nao=nao,nmo=nmo,&
|
||||
uniform_occupation=uni_occ,maxocc=maxocc)
|
||||
|
||||
NULLIFY(fm_work1)
|
||||
! create a new matrix
|
||||
NULLIFY(tmp_fm_struct)
|
||||
CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,&
|
||||
ncol_global=nmo,para_env=para_env,context=mo_coeff%matrix_struct%context,error=error)
|
||||
IF(ASSOCIATED(fm_work1))THEN
|
||||
CALL cp_fm_release(fm_work1,error=error)
|
||||
END IF
|
||||
CALL cp_fm_create (fm_work1, tmp_fm_struct ,error=error)
|
||||
CALL cp_fm_set_all(fm_work1,0.0_dp,error=error)
|
||||
IF ( .NOT. uni_occ ) THEN
|
||||
CALL cp_fm_create (psi0_fk, tmp_fm_struct ,error=error)
|
||||
CALL cp_fm_set_all(psi0_fk,0.0_dp,error=error)
|
||||
ELSE
|
||||
NULLIFY(psi0_fk)
|
||||
END IF
|
||||
CALL cp_fm_struct_release ( tmp_fm_struct ,error=error)
|
||||
CALL cp_fm_release(fm_work1,error=error)
|
||||
IF(.NOT.uni_occ) THEN
|
||||
CALL cp_fm_release(psi0_fk,error=error)
|
||||
ELSE
|
||||
NULLIFY(psi0_fk)
|
||||
ENDIF
|
||||
ENDDO ! ispin
|
||||
|
||||
! Allocate buffer vectors
|
||||
DO idir = 1,3
|
||||
ALLOCATE(vecbuf(idir)%array(1,nao),STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
END DO
|
||||
ALLOCATE(vecbuf_psi1(1,nao),STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
! First compute the contribution to the CHI tensor coming from the soft density
|
||||
CALL chi_soft_analytic(nmr_env,qs_env,ispin,i_B,error=error)
|
||||
|
||||
! Construct the 3 density matrices for the field in direction i_B
|
||||
|
||||
! First the full matrix psi_a_iB
|
||||
! psi1_rxp - psi1_D - (R_\nu-dk)_ii psi1_piiiB + (R_\nu-dk)_iii psi1_piiB
|
||||
psi_a_iB => psi1(ispin)%matrix
|
||||
CALL cp_fm_set_all(psi_a_iB,0.0_dp,error=error)
|
||||
|
||||
! contributions from the response psi1_p_ii and psi1_p_iii
|
||||
DO istate = 1,nmr_env%nstates(ispin)
|
||||
CALL cp_fm_get_submatrix(nmr_env%psi1_p(ispin,ii_B)%matrix,&
|
||||
vecbuf(ii_B)%array,1,istate,nao,1,transpose=.TRUE.,&
|
||||
error=error)
|
||||
CALL cp_fm_get_submatrix(nmr_env%psi1_p(ispin,iii_B)%matrix,&
|
||||
vecbuf(iii_B)%array,1,istate,nao,1,transpose=.TRUE.,&
|
||||
error=error)
|
||||
dk(1:3) = nmr_env%centers_set(ispin)%array(1:3,istate)
|
||||
DO iao = 1,nao
|
||||
Rmu(1:3) = nmr_env%basisfun_center(1:3,iao)
|
||||
Rmu_dk = pbc(dk,Rmu,cell)
|
||||
vecbuf_psi1(1,iao) = Rmu_dk(ii_B)*vecbuf(iii_B)%array(1,iao)-&
|
||||
Rmu_dk(iii_B)*vecbuf(ii_B)%array(1,iao)
|
||||
END DO ! iao
|
||||
! Copy the vector in the full matrix psi1
|
||||
CALL cp_fm_set_submatrix(psi1(ispin)%matrix,vecbuf_psi1,&
|
||||
1,istate,nao,1,transpose=.TRUE.,error=error)
|
||||
END DO ! istate
|
||||
|
||||
|
||||
! contribution from the response psi1_rxp
|
||||
CALL cp_fm_scale_and_add(-1.0_dp,psi_a_iB,1.0_dp,&
|
||||
nmr_env%psi1_rxp(ispin,i_B)%matrix,error=error)
|
||||
IF(nmr_env%full_nmr) THEN
|
||||
! contribution from the response psi1_D
|
||||
CALL cp_fm_scale_and_add(1.0_dp,psi_a_iB,-1.0_dp,&
|
||||
nmr_env%psi1_D(ispin,i_B)%matrix,error=error)
|
||||
END IF
|
||||
|
||||
! Multiply by the occupation number for the density matrix
|
||||
IF ( .NOT. uni_occ ) THEN
|
||||
alpha = 1.0_dp
|
||||
CALL cp_fm_to_fm(mo_coeff,psi0_fk,error=error)
|
||||
CALL cp_fm_column_scale(psi0_fk,occupation(1:homo))
|
||||
ELSE
|
||||
alpha = maxocc
|
||||
psi0_fk => mo_coeff
|
||||
END IF
|
||||
|
||||
! Build the first density matrix
|
||||
density_matrix(ispin)%matrix => nmr_env%jp1_ao
|
||||
CALL set_matrix(density_matrix(ispin)%matrix ,0.0_dp)
|
||||
CALL cp_sm_plus_fm_fm_t(sparse_matrix=density_matrix(ispin)%matrix ,&
|
||||
matrix_v=psi0_fk,&
|
||||
matrix_g=psi_a_iB,&
|
||||
ncol=homo,&
|
||||
alpha=alpha,error=error)
|
||||
|
||||
! Build the second density matrix
|
||||
density_matrix_iii(ispin)%matrix => nmr_env%jp2_ao(1)%matrix
|
||||
CALL set_matrix(density_matrix_iii(ispin)%matrix ,0.0_dp)
|
||||
CALL cp_sm_plus_fm_fm_t(sparse_matrix=density_matrix_iii(ispin)%matrix ,&
|
||||
matrix_v=psi0_fk,&
|
||||
matrix_g=nmr_env%psi1_p(ispin,iii_B)%matrix,&
|
||||
ncol=homo,&
|
||||
alpha=alpha,error=error)
|
||||
|
||||
! Build the third density matrix
|
||||
density_matrix_ii(ispin)%matrix => nmr_env%jp2_ao(2)%matrix
|
||||
CALL set_matrix(density_matrix_ii(ispin)%matrix ,0.0_dp)
|
||||
CALL cp_sm_plus_fm_fm_t(sparse_matrix=density_matrix_ii(ispin)%matrix ,&
|
||||
matrix_v=psi0_fk,&
|
||||
matrix_g=nmr_env%psi1_p(ispin,ii_B)%matrix,&
|
||||
ncol=homo,&
|
||||
alpha=alpha,error=error)
|
||||
|
||||
DO idir3=1,3
|
||||
! set to zero for the calculation of the shift
|
||||
CALL coeff_zero(shift_pw_gspace(idir3))
|
||||
END DO
|
||||
DO idir = 1,3
|
||||
|
||||
! Calculate the current density on the pw grid (only soft if GAPW)
|
||||
! idir is the cartesian component of the response current density
|
||||
! generated by the magnetic field pointing in cartesian direction i_B
|
||||
! Use the qs_rho_type already used for rho during the scf
|
||||
IF(gapw) THEN
|
||||
DO idir = 1,3
|
||||
! compute the atomic response current densities on the spherical grids
|
||||
! First the sparse matrices are multiplied by the expansion coefficients
|
||||
! this is the eqiuvalent of the CPC for the charge density
|
||||
CALL calculate_jrho_atom_coeff(qs_env,nmr_env,density_matrix,&
|
||||
& density_matrix_ii,density_matrix_iii,&
|
||||
& i_B,idir,error=error)
|
||||
! Then the radial parts are computed on the local radial grid, atom by atom
|
||||
! 8 functions are computed for each atom, per grid point
|
||||
! and per LM angular momentum. The multiplication by the Clebsh-Gordon
|
||||
! coefficients or they correspondent for the derivatives, is also done here
|
||||
CALL calculate_jrho_atom_rad(qs_env,nmr_env,I_B,idir,error=error)
|
||||
|
||||
IF(nmr_env%store_current) THEN
|
||||
rho_rspace => nmr_env%jrho1_set(idir,i_B)%rho%rho_r(ispin)
|
||||
rho_gspace => nmr_env%jrho1_set(idir,i_B)%rho%rho_g(ispin)
|
||||
ELSE
|
||||
rho_rspace => rho_struct%rho_r(ispin)
|
||||
rho_gspace => rho_struct%rho_g(ispin)
|
||||
END IF
|
||||
! Finally the radial functions are multiplied by the YLM and properly summed
|
||||
! The resulting array is J on the local grid. One array per atom.
|
||||
|
||||
CALL coeff_zero(rho_rspace)
|
||||
CALL coeff_zero(rho_gspace)
|
||||
CALL calculate_jrho_resp(density_matrix(ispin)%matrix , &
|
||||
density_matrix_ii(ispin)%matrix ,&
|
||||
density_matrix_iii(ispin)%matrix , &
|
||||
i_B, idir, rho_rspace, rho_gspace,qs_env,&
|
||||
gapw, error=error)
|
||||
! Local contributions by numerical integration over the spherical grids
|
||||
CALL shift_atom(qs_env,nmr_env,i_B,idir,error=error)
|
||||
ENDDO ! idir
|
||||
ENDIF
|
||||
!
|
||||
! Tranform the column i_B of the shift tensor from reciprocal to real space
|
||||
! get the values on the atomic positions and on the other required points (if any)
|
||||
! this can be done by interpolation of the values of the grid
|
||||
! on the required positions in real space
|
||||
DO idir = 1,3
|
||||
! Here the shifts are computed from summation of the coeff on the G-grip .
|
||||
CALL gsum_shift_pwgrid(nmr_env,pw_env,particle_set,cell,&
|
||||
& shift_pw_gspace(idir),i_B,idir,nmr_section,&
|
||||
& error=error)
|
||||
! Here first G->R and then interpolation to get the shifts.
|
||||
! The interpolation doesnt work in parallel yet.
|
||||
! The choice between both methods should be left to the user.
|
||||
!CALL coeff_transform_space(shift_pw_gspace(idir),shift_pw_rspace)
|
||||
!CALL interpolate_shift_pwgrid(nmr_env,pw_env,particle_set,cell,shift_pw_rspace,&
|
||||
! i_B,idir,nmr_section,error=error)
|
||||
ENDDO
|
||||
|
||||
scale_fac = cell%deth / twopi
|
||||
CALL coeff_scale(rho_rspace,scale_fac)
|
||||
CALL coeff_scale(rho_gspace,scale_fac)
|
||||
IF (BTEST(cp_print_key_should_output(logger%iter_info,nmr_section,&
|
||||
& "PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN
|
||||
|
||||
jrho_tot_G(idir,i_B) = calculate_total_rho(rho_gspace)
|
||||
jrho_tot_R(idir,i_B) = calculate_total_rho(rho_rspace)
|
||||
|
||||
IF(output_unit>0) THEN
|
||||
WRITE(output_unit,'(T2,2(A,E24.16))') 'Integrated j_'&
|
||||
&//ACHAR(idir+119)//ACHAR(i_B+119)//'(r): G-space=',&
|
||||
&jrho_tot_G(idir,i_B),' R-space=',jrho_tot_R(idir,i_B)
|
||||
END IF
|
||||
!
|
||||
!Field gradient
|
||||
! loop over the Gvec components: x,y,z
|
||||
DO idir2 = 1,3
|
||||
IF(idir /= idir2) THEN
|
||||
! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i))
|
||||
|
||||
my_chi = 0.0_dp !nmr_env%chi_SI2shiftppm/nmr_env%shift_factor*&
|
||||
! nmr_env%chi_tensor(idir3,i_B)
|
||||
CALL mult_G_ov_G2_grid(cell,auxbas_pw_pool,rho_gspace,&
|
||||
pw_gspace_work,idir2,my_chi,error=error)
|
||||
|
||||
! scale and add to the correct component of the shift column
|
||||
CALL set_vecp_rev(idir,idir2,idir3)
|
||||
scale_fac=fac_vecp(idir3,idir2,idir)
|
||||
CALL coeff_scale(pw_gspace_work,scale_fac)
|
||||
|
||||
CALL coeff_sumup(pw_gspace_work,shift_pw_gspace(idir3))
|
||||
END IF
|
||||
END DO
|
||||
END DO ! idir
|
||||
!
|
||||
! Deallocate buffer vectors
|
||||
DO idir = 1,3
|
||||
DEALLOCATE(vecbuf(idir)%array,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
END DO
|
||||
DEALLOCATE(vecbuf_psi1,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
CALL cp_fm_release(fm_work1,error=error)
|
||||
IF(.NOT.uni_occ) THEN
|
||||
CALL cp_fm_release(psi0_fk,error=error)
|
||||
ELSE
|
||||
NULLIFY(psi0_fk)
|
||||
ENDIF
|
||||
END DO ! ispin
|
||||
|
||||
IF(gapw) THEN
|
||||
DO idir = 1,3
|
||||
! compute the atomic response current densities on the spherical grids
|
||||
! First the sparse matrices are multiplied by the expansion coefficients
|
||||
! this is the eqiuvalent of the CPC for the charge density
|
||||
CALL calculate_jrho_atom_coeff(qs_env,nmr_env,&
|
||||
density_matrix, density_matrix_ii,&
|
||||
density_matrix_iii, i_B, idir, error=error)
|
||||
! Then the radial parts are computed on the local radial grid, atom by atom
|
||||
! 8 functions are computed for each atom, per grid point
|
||||
! and per LM angular momentum. The multiplication by the Clebsh-Gordon
|
||||
! coefficients or they correspondent for the derivatives, is also done here
|
||||
CALL calculate_jrho_atom_rad(qs_env,nmr_env,I_B,idir,error=error)
|
||||
|
||||
! Finally the radial functions are multiplied by the YLM and properly summed
|
||||
! The resulting array is J on the local grid. One array per atom.
|
||||
|
||||
! Local contributions by numerical integration over the spherical grids
|
||||
CALL shift_atom(qs_env,nmr_env,i_B,idir,error=error)
|
||||
END DO ! idir
|
||||
END IF
|
||||
|
||||
! Tranform the column i_B of the shift tensor from reciprocal to real space
|
||||
! get the values on the atomic positions and on the other required points (if any)
|
||||
! this can be done by interpolation of the values of the grid
|
||||
! on the required positions in real space
|
||||
|
||||
DO idir = 1,3
|
||||
CALL coeff_transform_space(shift_pw_gspace(idir),shift_pw_rspace)
|
||||
CALL interpolate_shift_pwgrid(nmr_env,pw_env,particle_set,cell,shift_pw_rspace,&
|
||||
i_B,idir,nmr_section,error=error)
|
||||
END DO
|
||||
|
||||
IF (BTEST(cp_print_key_should_output(logger%iter_info,nmr_section,&
|
||||
"PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN
|
||||
|
||||
DO idir = 1,3
|
||||
CALL coeff_zero(wf_r)
|
||||
!!WARNING
|
||||
! if nspins=2 here ist should be add instead of copy, or separated cubefiles
|
||||
! for the different spin should be written
|
||||
!!WARNING
|
||||
DO ispin =1 ,nspins
|
||||
CALL coeff_copy(nmr_env%jrho1_set(idir,i_B)%rho%rho_r(ispin),wf_r)
|
||||
END DO
|
||||
IF(gapw) THEN
|
||||
! Add the local hard and soft contributions
|
||||
! This can be done atom by atom by a spline extrapolation of the values
|
||||
! on the spherical grid to the grid points.
|
||||
DO iatom = 1,natom
|
||||
END DO
|
||||
END IF
|
||||
filename="jresp"
|
||||
WRITE(ext,'(a2,I1,a2,I1,a5)') "iB",i_B,"_d",idir,".cube"
|
||||
WRITE(ext,'(a2,a1,a2,a1,a5)') "iB",labels(i_B),"_d",labels(idir),".cube"
|
||||
unit_nr=cp_print_key_unit_nr(logger,nmr_section,"PRINT%CURRENT_CUBES",&
|
||||
extension=TRIM(ext),middle_name=TRIM(filename),&
|
||||
log_filename=.FALSE.,file_position="REWIND",error=error)
|
||||
|
||||
CALL rs_pw_to_cube(wf_r%pw,unit_nr,ionode,"RESPONSE CURRENT DENSITY ",&
|
||||
! rs,stride=section_get_ival(nmr_section,"PRINT%CURRENT_CUBES%STRIDE"),&
|
||||
stride=section_get_ival(nmr_section,"PRINT%CURRENT_CUBES%STRIDE",error=error),&
|
||||
error=error)
|
||||
CALL cp_print_key_finished_output(unit_nr,logger,nmr_section,&
|
||||
"PRINT%CURRENT_CUBES",error=error)
|
||||
END DO
|
||||
|
||||
END IF
|
||||
END DO ! i_B
|
||||
!
|
||||
! Integrated current response checksum
|
||||
IF(output_unit>0) THEN
|
||||
WRITE(output_unit,'(T2,A,E24.16)') 'CheckSum R-integrated j=',&
|
||||
SQRT(DDOT(9,jrho_tot_R(1,1),1,jrho_tot_R(1,1),1))
|
||||
ENDIF
|
||||
|
||||
IF (BTEST(cp_print_key_should_output(logger%iter_info,nmr_section,&
|
||||
"PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,wf_r,&
|
||||
error=error)
|
||||
CALL rs_pool_give_back_rs(auxbas_rs_pool,rs, error=error)
|
||||
END IF
|
||||
|
||||
!-----------------------------------------------------------------------!
|
||||
! Dellocate grids for the calculation of jrho and the shift
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,pw_gspace_work,&
|
||||
error=error)
|
||||
DO idir = 1,3
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,shift_pw_gspace(idir), &
|
||||
error=error)
|
||||
END DO
|
||||
DEALLOCATE(shift_pw_gspace,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,shift_pw_rspace,&
|
||||
error=error)
|
||||
!
|
||||
!-----------------------------------------------------------------------!
|
||||
DO ispin = 1,nspins
|
||||
NULLIFY(density_matrix(ispin)%matrix)
|
||||
NULLIFY(density_matrix_ii(ispin)%matrix)
|
||||
NULLIFY(density_matrix_iii(ispin)%matrix)
|
||||
END DO
|
||||
DEALLOCATE(density_matrix,density_matrix_ii,density_matrix_iii,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
|
||||
CALL timestop(0.0_dp,handle)
|
||||
END IF ! failure
|
||||
DO idir = 1,3
|
||||
CALL coeff_zero(wf_r)
|
||||
!!WARNING
|
||||
! if nspins=2 here ist should be add instead of copy, or separated cubefiles
|
||||
! for the different spin should be written
|
||||
!!WARNING
|
||||
DO ispin =1 ,nspins
|
||||
CALL coeff_copy(nmr_env%jrho1_set(idir,i_B)%rho%rho_r(ispin),wf_r)
|
||||
ENDDO
|
||||
IF(gapw) THEN
|
||||
! Add the local hard and soft contributions
|
||||
! This can be done atom by atom by a spline extrapolation of the values
|
||||
! on the spherical grid to the grid points.
|
||||
CALL stop_program(routineP,"GAPW needs to be finalized")
|
||||
DO iatom = 1,natom
|
||||
ENDDO
|
||||
ENDIF
|
||||
filename="jresp"
|
||||
WRITE(ext,'(a2,I1,a2,I1,a5)') "iB",i_B,"_d",idir,".cube"
|
||||
WRITE(ext,'(a2,a1,a2,a1,a5)') "iB",ACHAR(i_B+119),"_d",ACHAR(idir+119),".cube"
|
||||
unit_nr=cp_print_key_unit_nr(logger,nmr_section,"PRINT%CURRENT_CUBES",&
|
||||
& extension=TRIM(ext),middle_name=TRIM(filename),&
|
||||
& log_filename=.FALSE.,file_position="REWIND",&
|
||||
& error=error)
|
||||
|
||||
CALL rs_pw_to_cube(wf_r%pw,unit_nr,ionode,"RESPONSE CURRENT DENSITY ",&
|
||||
!rs,stride=section_get_ival(nmr_section,"PRINT%CURRENT_CUBES%STRIDE"),&
|
||||
& stride=section_get_ival(nmr_section,"PRINT%CURRENT_CUBES%STRIDE",&
|
||||
& error=error),error=error)
|
||||
CALL cp_print_key_finished_output(unit_nr,logger,nmr_section,&
|
||||
& "PRINT%CURRENT_CUBES",error=error)
|
||||
ENDDO
|
||||
|
||||
ENDIF ! current cube
|
||||
ENDDO ! i_B
|
||||
!
|
||||
! Integrated current response checksum
|
||||
IF(output_unit>0) THEN
|
||||
WRITE(output_unit,'(T2,A,E24.16)') 'CheckSum R-integrated j=',&
|
||||
SQRT(DDOT(9,jrho_tot_R(1,1),1,jrho_tot_R(1,1),1))
|
||||
ENDIF
|
||||
!
|
||||
IF (BTEST(cp_print_key_should_output(logger%iter_info,nmr_section,&
|
||||
& "PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,wf_r,error=error)
|
||||
CALL rs_pool_give_back_rs(auxbas_rs_pool,rs, error=error)
|
||||
ENDIF
|
||||
!
|
||||
! Dellocate grids for the calculation of jrho and the shift
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,pw_gspace_work,error=error)
|
||||
DO idir = 1,3
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,shift_pw_gspace(idir),error=error)
|
||||
ENDDO
|
||||
DEALLOCATE(shift_pw_gspace,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
CALL pw_pool_give_back_coeff(auxbas_pw_pool,shift_pw_rspace,error=error)
|
||||
!
|
||||
!
|
||||
DO ispin = 1,nspins
|
||||
NULLIFY(density_matrix(ispin)%matrix)
|
||||
NULLIFY(density_matrix_ii(ispin)%matrix)
|
||||
NULLIFY(density_matrix_iii(ispin)%matrix)
|
||||
ENDDO
|
||||
DEALLOCATE(density_matrix,density_matrix_ii,density_matrix_iii,STAT=istat)
|
||||
CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
|
||||
ENDIF ! failure
|
||||
!
|
||||
! Finalize
|
||||
CALL timestop(0.0_dp,handle)
|
||||
!
|
||||
END SUBROUTINE nmr_response_current
|
||||
|
||||
! *****************************************************************************
|
||||
|
|
|
|||
|
|
@ -122,7 +122,7 @@ CONTAINS
|
|||
nsgf, output_unit
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
|
||||
LOGICAL :: failure, ionode
|
||||
REAL(dp) :: ck(3), ckdk(3), dk(3)
|
||||
REAL(dp) :: ck(3), ckdk(3), dk(3), chk(3)
|
||||
REAL(dp), DIMENSION(:, :), POINTER :: basisfun_center, vecbuf_c0
|
||||
TYPE(atomic_kind_type), DIMENSION(:), &
|
||||
POINTER :: atomic_kind_set
|
||||
|
|
@ -147,10 +147,6 @@ CONTAINS
|
|||
POINTER :: op_p_ao, op_rmd_ao
|
||||
TYPE(section_vals_type), POINTER :: lr_section
|
||||
|
||||
real*8,allocatable,dimension(:,:)::buffer1,buffer2,buffer3
|
||||
!dbg
|
||||
!dbg
|
||||
|
||||
CALL timeset(routineN,"I"," ",handle)
|
||||
|
||||
failure = .FALSE.
|
||||
|
|
@ -201,9 +197,11 @@ CONTAINS
|
|||
!!$ !First term: operator matrix elements
|
||||
!!$ CALL rmc_x_p_xyz_ao(op_rmd_ao,qs_env,minimum_image=.FALSE.,error=error)
|
||||
!************************************************************
|
||||
write(*,*) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING'
|
||||
write(*,*) 'rmc_x_p_xyz_ao is computed in a different way'
|
||||
write(*,*) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING'
|
||||
IF(output_unit>0) THEN
|
||||
WRITE(output_unit,*) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING'
|
||||
WRITE(output_unit,*) 'rmc_x_p_xyz_ao: L is computed for each states'
|
||||
WRITE(output_unit,*) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING'
|
||||
ENDIF
|
||||
|
||||
DO ispin=1,dft_control%nspins
|
||||
p_psi0(ispin,1)%matrix%local_data=0d0
|
||||
|
|
@ -232,27 +230,27 @@ CONTAINS
|
|||
ENDDO
|
||||
!
|
||||
! print checksums
|
||||
chk(1) = checksum_matrix(op_rmd_ao(1)%matrix,para_env)
|
||||
chk(2) = checksum_matrix(op_rmd_ao(2)%matrix,para_env)
|
||||
chk(3) = checksum_matrix(op_rmd_ao(3)%matrix,para_env)
|
||||
IF(output_unit>0) THEN
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum L_x =', &
|
||||
& checksum_matrix(op_rmd_ao(1)%matrix)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum L_y =', &
|
||||
& checksum_matrix(op_rmd_ao(2)%matrix)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum L_z =', &
|
||||
& checksum_matrix(op_rmd_ao(3)%matrix)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum L_x =',chk(1)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum L_y =',chk(2)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum L_z =',chk(3)
|
||||
ENDIF
|
||||
!************************************************************
|
||||
!
|
||||
! Calculate the px py pz operators
|
||||
CALL p_xyz_ao(op_p_ao,qs_env,minimum_image=.FALSE.,error=error)
|
||||
|
||||
!
|
||||
! print checksums
|
||||
chk(1) = checksum_matrix(op_p_ao(1)%matrix,para_env)
|
||||
chk(2) = checksum_matrix(op_p_ao(2)%matrix,para_env)
|
||||
chk(3) = checksum_matrix(op_p_ao(3)%matrix,para_env)
|
||||
IF(output_unit>0) THEN
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum P_x =', &
|
||||
& checksum_matrix(op_p_ao(1)%matrix)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum P_y =', &
|
||||
& checksum_matrix(op_p_ao(2)%matrix)
|
||||
write(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum P_z =', &
|
||||
& checksum_matrix(op_p_ao(3)%matrix)
|
||||
WRITE(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum P_x =',chk(1)
|
||||
WRITE(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum P_y =',chk(2)
|
||||
WRITE(output_unit,'(A,E23.16)') 'nmr_operators: CheckSum P_z =',chk(3)
|
||||
ENDIF
|
||||
! Apply the p operator to the psi0
|
||||
DO idir = 1,3
|
||||
|
|
|
|||
|
|
@ -94,18 +94,20 @@ MODULE qs_linres_nmr_shift
|
|||
real_matrix_p_type,&
|
||||
replicate_matrix_structure,&
|
||||
set_matrix
|
||||
USE message_passing, ONLY: mp_sum
|
||||
USE termination, ONLY: stop_program
|
||||
USE timings, ONLY: timeset,&
|
||||
timestop
|
||||
|
||||
#include "cp_common_uses.h"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
|
||||
PRIVATE
|
||||
|
||||
! *** Public subroutines ***
|
||||
PUBLIC :: chi_soft_analytic, interpolate_shift_pwgrid, mult_G_ov_G2_grid, nmr_print_shift
|
||||
PUBLIC :: chi_soft_analytic, interpolate_shift_pwgrid, &
|
||||
& gsum_shift_pwgrid, mult_G_ov_G2_grid, nmr_print_shift
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_nmr_shift'
|
||||
|
||||
|
|
@ -151,12 +153,9 @@ CONTAINS
|
|||
TYPE(real_matrix_p_type), &
|
||||
DIMENSION(:, :), POINTER :: op_mom_der_ao
|
||||
|
||||
!dbg
|
||||
! TYPE(cp_fm_type), POINTER :: fm_work2
|
||||
!dbg
|
||||
|
||||
CALL timeset(routineN,"I"," ",handle)
|
||||
|
||||
|
||||
failure = .FALSE.
|
||||
|
||||
NULLIFY(dft_control,matrix_s,mos,para_env)
|
||||
|
|
@ -483,12 +482,12 @@ CONTAINS
|
|||
!!*** **********************************************************************
|
||||
|
||||
SUBROUTINE interpolate_shift_pwgrid(nmr_env,pw_env,particle_set,cell,shift_pw_rspace,&
|
||||
i_B,idir,nmr_section,error)
|
||||
& i_B,idir,nmr_section,error)
|
||||
|
||||
TYPE(nmr_env_type) :: nmr_env
|
||||
TYPE(pw_env_type), POINTER :: pw_env
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particle_set
|
||||
POINTER :: particle_set
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
TYPE(coeff_type) :: shift_pw_rspace
|
||||
INTEGER, INTENT(IN) :: i_B, idir
|
||||
|
|
@ -496,8 +495,7 @@ CONTAINS
|
|||
TYPE(cp_error_type), INTENT(INOUT) :: error
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'interpolate_shift_pwgrid', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
routineP = moduleN//':'//routineN
|
||||
INTEGER :: aint_precond, iat, iatom, &
|
||||
max_iter, n_nics, natom, &
|
||||
precond_kind
|
||||
|
|
@ -513,8 +511,9 @@ CONTAINS
|
|||
TYPE(pw_spline_precond_type), POINTER :: precond
|
||||
TYPE(section_vals_type), POINTER :: interp_section
|
||||
|
||||
!dbg
|
||||
!dbg
|
||||
#ifdef __parallel
|
||||
CALL stop_program(routineP,"Doesnt support parallel yet.")
|
||||
#endif
|
||||
|
||||
failure = .FALSE.
|
||||
NULLIFY (interp_section)
|
||||
|
|
@ -563,11 +562,11 @@ CONTAINS
|
|||
ENDIF
|
||||
|
||||
DO iat = 1,natom
|
||||
iatom = cs_atom_list(iat)
|
||||
R_iatom = pbc(particle_set(iatom)%r,cell)
|
||||
shift_val = Eval_Interp_Spl3_pbc(R_iatom,shiftspl%pw,error=error)
|
||||
chemical_shift(idir,i_B,iatom)= chemical_shift(idir,i_B,iatom)+&
|
||||
nmr_env%shift_factor*shift_val
|
||||
iatom = cs_atom_list(iat)
|
||||
R_iatom = pbc(particle_set(iatom)%r,cell)
|
||||
shift_val = Eval_Interp_Spl3_pbc(R_iatom,shiftspl%pw,error=error)
|
||||
chemical_shift(idir,i_B,iatom)= chemical_shift(idir,i_B,iatom)+&
|
||||
nmr_env%shift_factor*shift_val
|
||||
END DO
|
||||
|
||||
IF(nmr_env%do_nics) THEN
|
||||
|
|
@ -586,6 +585,81 @@ CONTAINS
|
|||
|
||||
END SUBROUTINE interpolate_shift_pwgrid
|
||||
|
||||
SUBROUTINE gsum_shift_pwgrid(nmr_env,pw_env,particle_set,cell,shift_pw_gspace,&
|
||||
& i_B,idir,nmr_section,error)
|
||||
TYPE(nmr_env_type) :: nmr_env
|
||||
TYPE(pw_env_type), POINTER :: pw_env
|
||||
TYPE(particle_type), DIMENSION(:), &
|
||||
POINTER :: particle_set
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
TYPE(coeff_type) :: shift_pw_gspace
|
||||
INTEGER, INTENT(IN) :: i_B, idir
|
||||
TYPE(section_vals_type), POINTER :: nmr_section
|
||||
TYPE(cp_error_type), INTENT(INOUT) :: error
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'gsum_shift_pwgrid', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: iat, iatom, n_nics, natom
|
||||
INTEGER, DIMENSION(:), POINTER :: cs_atom_list
|
||||
LOGICAL :: failure
|
||||
REAL(dp) :: R_iatom(3), &
|
||||
ra(3)
|
||||
REAL(dp), DIMENSION(:, :), POINTER :: r_nics
|
||||
REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, &
|
||||
chemical_shift_nics
|
||||
!
|
||||
failure = .FALSE.
|
||||
NULLIFY(cs_atom_list,chemical_shift,chemical_shift_nics,r_nics)
|
||||
CPPrecondition(ASSOCIATED(shift_pw_gspace%pw),cp_failure_level,routineP,error,failure)
|
||||
IF(.NOT. failure) THEN
|
||||
!
|
||||
CALL get_nmr_env(nmr_env=nmr_env, cs_atom_list=cs_atom_list,&
|
||||
& chemical_shift=chemical_shift,&
|
||||
& chemical_shift_nics=chemical_shift_nics,&
|
||||
& n_nics=n_nics,r_nics=r_nics,error=error)
|
||||
!
|
||||
IF(ASSOCIATED(cs_atom_list)) THEN
|
||||
natom = SIZE(cs_atom_list,1)
|
||||
ELSE
|
||||
natom = -1
|
||||
ENDIF
|
||||
!
|
||||
! compute the chemical shift
|
||||
DO iat = 1,natom
|
||||
iatom = cs_atom_list(iat)
|
||||
R_iatom = pbc(particle_set(iatom)%r,cell)
|
||||
chemical_shift(idir,i_B,iatom)= chemical_shift(idir,i_B,iatom)+&
|
||||
& nmr_env%shift_factor*REAL(gsumr(R_iatom,shift_pw_gspace%pw),dp)
|
||||
ENDDO
|
||||
!
|
||||
! compute nics
|
||||
IF(nmr_env%do_nics) THEN
|
||||
DO iat = 1,n_nics
|
||||
ra = pbc(r_nics(:,iat),cell)
|
||||
chemical_shift_nics(idir,i_B,iatom)= chemical_shift_nics(idir,i_B,iatom)+&
|
||||
& nmr_env%shift_factor*REAL(gsumr(ra,shift_pw_gspace%pw),dp)
|
||||
ENDDO
|
||||
ENDIF
|
||||
ENDIF
|
||||
END SUBROUTINE gsum_shift_pwgrid
|
||||
|
||||
FUNCTION gsumr(r,pw)
|
||||
REAL(dp), INTENT(IN) :: r(3)
|
||||
TYPE(pw_type), POINTER :: pw
|
||||
TYPE(pw_grid_type), POINTER :: grid
|
||||
COMPLEX(dp) :: gsumr,rg
|
||||
INTEGER :: ig
|
||||
grid => pw%pw_grid
|
||||
gsumr = CMPLX(0.0_dp,0.0_dp,KIND=dp)
|
||||
DO ig = grid%first_gne0,grid%ngpts_cut_local
|
||||
rg = ( grid%g(1,ig)*r(1) + grid%g(2,ig)*r(2) + grid%g(3,ig)*r(3) ) * gaussi
|
||||
gsumr = gsumr + pw%cc(ig)*EXP(rg)
|
||||
ENDDO
|
||||
IF(grid%have_g0) gsumr = gsumr + pw%cc(1)
|
||||
CALL mp_sum(gsumr,grid%para%group)
|
||||
END FUNCTION gsumr
|
||||
|
||||
!!****f* qs_linres_nmr_shift%mult_G_ov_G2_grid
|
||||
!!
|
||||
!! NAME
|
||||
|
|
|
|||
|
|
@ -2117,8 +2117,9 @@ END SUBROUTINE cp_sm_scale_and_add
|
|||
|
||||
! *****************************************************************************
|
||||
|
||||
FUNCTION checksum_real_matrix(a)
|
||||
FUNCTION checksum_real_matrix(a,para_env)
|
||||
TYPE(real_matrix_type), POINTER :: a
|
||||
TYPE(cp_para_env_type), POINTER :: para_env
|
||||
CHARACTER(LEN=*), PARAMETER :: &
|
||||
routine = "SUBROUTINE checksum_real_matrix (MODULE sparse_matrix_types)"
|
||||
INTEGER :: iblock_row,n
|
||||
|
|
@ -2132,6 +2133,9 @@ END SUBROUTINE cp_sm_scale_and_add
|
|||
IF (.NOT.ASSOCIATED(a)) THEN
|
||||
CALL stop_program(routine,"The matrix pointer a is not associated")
|
||||
END IF
|
||||
IF (.NOT.ASSOCIATED(para_env)) THEN
|
||||
CALL stop_program(routine,"The para_env pointer a is not associated")
|
||||
END IF
|
||||
DO iblock_row=1,a%nblock_row
|
||||
block_node => first_block_node(a,iblock_row)
|
||||
DO WHILE (ASSOCIATED(block_node))
|
||||
|
|
@ -2142,6 +2146,7 @@ END SUBROUTINE cp_sm_scale_and_add
|
|||
block_node => block_node%next_block_node
|
||||
END DO
|
||||
END DO
|
||||
CALL mp_sum(checksum_real_matrix,para_env%group)
|
||||
END FUNCTION checksum_real_matrix
|
||||
! *****************************************************************************
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue