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

[Laino Teodoro]

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