[CP2K:1932] Re: QMMM: MM SPLINE: no convergence
Peter Mamonov
pmam... at gmail.com
Mon Mar 23 11:39:37 UTC 2009
Thank you for answers, Teodoro : )
Peter
On Mon, Mar 23, 2009 at 2:51 AM, Laino Teodoro <teodor... at gmail.com> wrote:
>
> Mechanical coupling means that the QM and the MM system interact
> electrostatically only at the classical level.
> If you want an electrostatic-coupling between QM and MM (where the
> effect of the charge is directly in the DFT scheme)
> you need to use GAUSSIAN for E_COUPL keyword.
>
> I would suggest you to read one of the main feature articles on QM/MM
> (about the possible coupling schemes).
>
> If you decide to go for the E_COUPL GAUSSIAN I would suggest you to
> read as well the paper on which
> it is based the implementation in CP2K.
>
> On 22 Mar 2009, at 11:55, Peter Mamonov wrote:
>
>>
>> Thank you Teodoro!
>>
>> Excuse me for bothering you, but I'm still in doubt. For example about
>> this piece of output:
>>
>> SPLINE_INFO| Spline number: 1
>> SPLINE_INFO| Spline number 1 Unique number 102 Kinds
>> involved: 1 1
>> SPLINE_INFO| Number of points used in the Splines: 545789
>> SPLINE_INFO| Accuracy requested: 0.000000100000
>> SPLINE_INFO| Accuracy achieved: 0.000000092134
>> SPLINE_INFO| RMIN and RMAX used for spline (bohr):
>> 0.900000 1000.000000
>> SPLINE_INFO| RMIN and RMAX used to achieve spline accuracy (bohr):
>> 0.900366 1000.000000
>> SPLINE_INFO| Spline value at RMIN (Hartree): 0.002769639
>> SPLINE_INFO| Spline value at RMAX (Hartree): 0.000000000
>> SPLINE_INFO| Electrostatic and Non-Bonded Cutoffs: 0.000002500
>> 0.000000000
>>
>> The code uses the value of ~500A (1e3 Bohr) for RMAX, but i've set 50A
>> cutoff for non-bonded interactions, and the system is non-periodic. It
>> also uses enormous number of points to define spline, compared to
>> other pairs.
>>
>> Do these figures look reasonable to you?
>
> This happens for the following reason:
>
> The van der waals between the QM types are switched off (the pair
> potential is then zeroed and by default the maximum upper limit is
> chosen, which is 1000.).
> Normally for electrostatic embedding this is not a problem since on
> the interval 0-1000.
> the spline is uniformly zero, so if you choose 1000. or 50. you will
> always use 2 points to
> map that spline.
> But if you choose (like you did) electrostatic coupling NONE
> (mechanical) then
> on this spline we need still to map the pure coulomb (1/r).
> While for pure FIST this potential is somehow attenuated by the LJ
> (therefore you see
> a better spline build-up), for QMMM runs the single Coulomb potential
> (1/r) is much difficult
> to converge not for the large cutoff (1000.) but for the function
> behavior at r->0.
> Keep in mind, in fact, that we don't use equi-spaced splines in the X-
> domain, but splines equispaced in the
> 1/X-domain (which means that roughly 90% of the points are
> accumulated for small values of R, close to the
> divergence of the 1/r potential).
>
> After all this.. I must say that what you observe is the correct
> behavior.
>
> Although playing with EPS* and other keywords may solve this issue..
> you should think more whether few things like
> the periodicity NONE for the classical part and e-coupling for qm/mm
> are really the things you want..
>
> Cheers,
> Teo
>
>>
>>
>>> Keep in mind that you are using the default E_COUPL, which is the
>>> mechanical coupling (the electrostatic QM-MM
>>> coupling is done at the level of classical forcefield).
>>
>> So, do i understand correctly: this means, that electrostatic term due
>> to MM point charges is added to the DFT functional, while MM atoms
>> interact with QM-generated point charges on QM atoms?
>>
>>
>> Thanks,
>> Peter
>>
>> On Sun, Mar 22, 2009 at 9:34 AM, Laino Teodoro
>> <teodor... at gmail.com> wrote:
>>>
>>> Dear Peter,
>>>
>>> what you see is the right behavior. On the spline it is mapped part
>>> of the FF (this means LJ or other similar
>>> potential + real-space part of coulomb).
>>> Since some of the above interactions are zeroed in the FF for QM-MM
>>> runs it is highly possible that some splines
>>> may have different behavior (different boundaries or different
>>> convergence).
>>>
>>> Keep in mind that you are using the default E_COUPL, which is the
>>> mechanical coupling (the electrostatic QM-MM
>>> coupling is done at the level of classical forcefield).
>>>
>>> Cheers,
>>> Teo
>>>
>>> On 21 Mar 2009, at 21:38, Peter wrote:
>>>
>>>>
>>>> Dear all,
>>>>
>>>> I'm trying to setup qm/mm calculations of a test system. MM
>>>> treatment
>>>> of a system using FIST works fine. However, QM/MM calculation fails
>>>> with the following message:
>>>>
>>>> *** 22:36:04 ERRORL2 in pair_potential:generate_spline_low
>>>> processor ***
>>>> *** 0 err=-300 SPLINE_INFO| Number of points: 2346797
>>>> obtained ***
>>>> *** accuracy 0.195114E-06. MM SPLINE: no convergence on required
>>>> accuracy ***
>>>> *** (adjust EPS_SPLINE and rerun) pair_potential.F line
>>>> 551 ***
>>>>
>>>> Adjusting EPS_SPLINE and/or R0_NB, EMAX_SPLINE solves the problem.
>>>> However I found this situation quite strange, since MM calculations
>>>> using the same FF and coordinates don't need such adjustments. Also
>>>> I've noticed, that during QMMM run spline generator uses much wider
>>>> boundaries (RMIN, RMAX) for some pairs, compared to MM run.
>>>>
>>>> So, is it the supposed behavior of the program? And if so, could you
>>>> comment on it, please?
>>>>
>>>> Thanks,
>>>> Peter
>>>>
>>>> Input and log files: http://groups.google.com/group/cp2k/web/
>>>> mm_spline-no_convergence.tar.gz
>>>>
>>>> # u10-qmmm.inp:
>>>> &GLOBAL
>>>> PROJECT U10-QMMM_OPT
>>>> RUN_TYPE GEO_OPT
>>>> PREFERRED_FFT_LIBRARY FFTMKL
>>>> &END GLOBAL
>>>>
>>>> &MOTION
>>>> &GEO_OPT
>>>> TYPE MINIMIZATION
>>>> MAX_ITER 10000
>>>> OPTIMIZER LBFGS
>>>> MAX_FORCE 5.0E-4
>>>> RMS_FORCE 2.0E-4
>>>> &LBFGS
>>>> MAX_H_RANK 10
>>>> MAX_F_PER_ITER 100
>>>> &END LBFGS
>>>> &END GEO_OPT
>>>>
>>>> &PRINT
>>>> &RESTART OFF
>>>> &END RESTART
>>>> &RESTART_HISTORY OFF
>>>> &END RESTART_HISTORY
>>>> &TRAJECTORY
>>>> FORMAT XYZ
>>>> FILENAME =u10-qmmm.xyz
>>>> &END TRAJECTORY
>>>> &END PRINT
>>>> &END MOTION
>>>>
>>>> &FORCE_EVAL
>>>> METHOD QMMM
>>>>
>>>> &DFT
>>>> CHARGE 0
>>>> MULTIPLICITY 1
>>>> BASIS_SET_FILE_NAME EMSL_BASIS_SETS
>>>> POTENTIAL_FILE_NAME POTENTIAL
>>>>
>>>> &MGRID
>>>> CUTOFF 250
>>>> COMMENSURATE
>>>> &END MGRID
>>>>
>>>> &QS
>>>> METHOD GAPW
>>>> &END QS
>>>>
>>>> &SCF
>>>> SCF_GUESS ATOMIC
>>>> EPS_SCF 1.0E-6
>>>> MAX_SCF 150
>>>> &END SCF
>>>>
>>>> &XC
>>>> #These are the coefficients used for B3LYP using VWN5, this is
>>>> recommended but doesn't match the default Gaussian definition
>>>> &XC_FUNCTIONAL
>>>> &LYP
>>>> SCALE_C 0.81
>>>> &END
>>>> &BECKE88
>>>> SCALE_X 0.72
>>>> &END
>>>> &VWN
>>>> FUNCTIONAL_TYPE VWN5
>>>> SCALE_C 0.19
>>>> &END
>>>> &XALPHA
>>>> SCALE_X 0.08
>>>> &END
>>>> &END XC_FUNCTIONAL
>>>> &HF
>>>> &SCREENING
>>>> EPS_SCHWARZ 1.0E-10
>>>> &END
>>>> &MEMORY
>>>> MAX_MEMORY 5
>>>> &END
>>>> FRACTION 0.20
>>>> &END
>>>> &END XC
>>>>
>>>> &POISSON
>>>> PERIODIC NONE
>>>> POISSON_SOLVER WAVELET
>>>> &END POISSON
>>>> &END DFT
>>>>
>>>>
>>>> &MM
>>>> &FORCEFIELD
>>>> PARM_FILE_NAME par_all27_prot_lipid.prm
>>>> PARMTYPE CHM
>>>> &SPLINE
>>>> RCUT_NB 50.0
>>>> # UNIQUE_SPLINE T
>>>> # EPS_SPLINE 1.0E-12
>>>> # EMAX_SPLINE 0.02
>>>> # R0_NB = 0.45
>>>> &END SPLINE
>>>> &END FORCEFIELD
>>>>
>>>> &POISSON
>>>> PERIODIC NONE
>>>> &EWALD
>>>> EWALD_TYPE NONE
>>>> &END EWALD
>>>> &END POISSON
>>>>
>>>> &NEIGHBOR_LISTS
>>>> GEO_CHECK F
>>>> &END NEIGHBOR_LISTS
>>>>
>>>> &PRINT
>>>> &FF_INFO
>>>> &END FF_INFO
>>>> &END PRINT
>>>> &END MM
>>>>
>>>> &QMMM
>>>> &CELL
>>>> ABC 15.0 15.0 15.0
>>>> &END CELL
>>>>
>>>> &QM_KIND H
>>>> MM_INDEX 14 15 16 17 18 19 20 21 22 28 29
>>>> &END QM_KIND
>>>>
>>>> &QM_KIND C
>>>> MM_INDEX 1 2 3 4 5 6 7 8 9 23
>>>> &END QM_KIND
>>>>
>>>> &QM_KIND O
>>>> MM_INDEX 10 11 12 13
>>>> &END QM_KIND
>>>>
>>>> &LINK
>>>> QM_INDEX 23
>>>> MM_INDEX 24
>>>> QM_KIND H
>>>> &END LINK
>>>>
>>>> &END QMMM
>>>>
>>>> &SUBSYS
>>>> &CELL
>>>> ABC 50.0 50.0 50.0
>>>> PERIODIC NONE
>>>> &END CELL
>>>>
>>>> &TOPOLOGY
>>>> COORD_FILE_NAME u10.pdb
>>>> COORD_FILE_FORMAT PDB
>>>> CONN_FILE_NAME u10.psf
>>>> CONN_FILE_FORMAT PSF
>>>> CENTER_COORDINATES
>>>> &END TOPOLOGY
>>>>
>>>> &KIND O
>>>> BASIS_SET 6-31G*
>>>> POTENTIAL ALL
>>>> &END KIND
>>>>
>>>> &KIND C
>>>> BASIS_SET 6-31G*
>>>> POTENTIAL ALL
>>>> &END KIND
>>>>
>>>> &KIND H
>>>> BASIS_SET 6-31G*
>>>> POTENTIAL ALL
>>>> &END KIND
>>>> &END SUBSYS
>>>> &END FORCE_EVAL
>>>>>
>>>
>>>
>>>>
>>>
>>
>> >
>
>
> >
>
More information about the CP2K-user
mailing list