Dear Amber Users,
We are studying the SCC-DFTB codes in AMBER 10, and had difficulties in understanding part of the code (we tried to understand the program based on the paper : J. Phys. Chem. A 111(26), 2007. We wonder if anyone could kindly help us out.
The first part that we do not understand is (qm2_dftb_scf.f, lines 481-503)
! COULOMBIC ELECTRONIC REPULSION ENERGY
! -------------------------------------
do i = 1, qmmm_struct%nquant_nlink
ecoul = ecoul + &
ks_struct%shift(i) &
*(qmat(i)+mcharge%qzero( izp_str%izp(i)))
end do
! EXTERNAL CHARGES (QM/MM) COULOMB ENERGY
! ---------------------------------------
! This term is a simple coulomb interaction between
! the external charge (in ks_struct%shiftE) and the mulliken charge
! on the atom.
! do i= qmmm_mpi%nquant_nlink_start, qmmm_mpi%nquant_nlink_end
do i= 1, qmmm_struct%nquant_nlink
eext = eext + ks_struct%shiftE(i)*(mcharge%qzero( izp_str%izp(i))-qmat(i))
end do
! Electronic Energy
! -----------------
! remark: elec_eng contains ks_struct%shiftE aready via ev,
! ks_struct%shift also contains -ks_struct%shiftE, i.e. ecoul also
! contains contributions from EXT
elec_eng = elec_eng-0.5d0*ecoul + 0.5d0*eext
My question for the code given above: why is this Ecoul subtracted from the electronic energy? Which equation does it correspond it to in the JPCA paper?
The second puzzle we had is lines 625 to 658 in qm2_dftb_scf.f,
! EWALD ENERGY CORRECTION
! -----------------------
! This correction is necessary because, up to this point,
! the energy for the Ewald sum only included half of the term from the MM atoms
! and half from the QM atoms. The QM atoms should contribute only half but the
! MM atoms should contribute full. The corrects qm_ewald_correct_ee routine
! supposedly corrects for this.
!
! That routine uses the density matrix, and expects it to be indexed according
! to the contents of qm2_params%orb_loc.
#ifndef SQM
if ( qmmm_nml%qm_ewald>0 ) then
! Calculate the correction in Hartree/Bohrs
call qm2_dftb_ewald_corr(qmmm_struct%nquant_nlink, ew_corr, &
qmewald%mmpot, scf_mchg)!qmat)
! Put correction into elec_eng.
elec_eng = elec_eng + ew_corr
end if
if ( qmmm_nml%qmgb == 2) then
!Calculate the correction for GB - removes double counting of
!GB energy in both escf and egb.
gb_escf_corr = 0.0d0
do i = 1,qmmm_struct%nquant_nlink
gb_escf_corr = gb_escf_corr + (qm_gb%gb_mmpot(i)+qm_gb%gb_qmpot(i))*scf_mchg(i)
end do
elec_eng = elec_eng + (0.5d0*gb_escf_corr*BOHRS_TO_A)
end if
#endif
Our question for this part is: why does QM part contribute only half and MM part contribute fully in the Ewald term? Will it affect the SCF calculations since it is now corrected out of the SCF loop?
Thank you very much for your help.
Best wishes,
Lijiang Yang
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Sep 20 2011 - 23:00:03 PDT
