[CP2K:1932] Re: QMMM: MM SPLINE: no convergence

[Peter Mamonov]

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
>>>>>
>>>
>>>
>>>>
>>>
>>
>> >
>
>
> >
>