Hi Lijiang,
Since I'm not very familiar with DFTB code, I'll just try to answer
your 2nd question.
To get QM/MM ewald energy, you would need two part of energy, QM-QM
energy and QM-MM energy. If we loop over QM-QM atom pairs, we'll count
the QM-QM energy twice. In contrast, if we loop over QM-MM pairs,
there'll be no double counting in the QM-MM energy. Please correct me
if I'm wrong here.
The details about the QM/MM ewald implementation can be found in the
paper: Ross et al. (2007) JCC, vol 29, p1019
Best,
Qiantao Wang
On Wed, Sep 21, 2011 at 1:45 AM, 杨立江 <lijiangy.pku.edu.cn> wrote:
>
> 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
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Sep 21 2011 - 15:00:03 PDT
