SCF, MD run-time verses atomic species.

350 views
Skip to first unread message

Simiam Ghan

unread,
Jul 1, 2016, 7:38:47 AM7/1/16
to cp2k
Dear all,
I am running NVT MD with Quickstep on a box of water with 64 molecules.   If I replace a water molecule with a KCl ion pair, i observe that the MD and SCF step times more than double on my setup.  (MD from ~20 sec to ~46 sec). SCF iterations are converging (except the very first MD step) in around 30 steps in both cases but each SCF now takes over twice as long as before.   Is there an explanation of why such a 'small' change in system could double the run time?  Is KCl really so heavy to calculate compared to H2O?  The number of electrons in the system increases by 8, from 512 to 520.  Below my KCl input file. 

Greetings,
Simiam



&GLOBAL

  PROJECT H2O_KCl
  RUN_TYPE MD
  PRINT_LEVEL MEDIUM

&END GLOBAL

&FORCE_EVAL

  METHOD Quickstep ! GPW method.

  &SUBSYS                                       ! A subsystem: coordinates, topology, molecules and cell.

    &CELL                                       ! Supercell setup.
      ABC [angstrom] 12.414 12.414 12.414 ! Using 64 H2O molecules, we thus get a density of 1g/cm^3.
      PERIODIC XYZ           ! Use PBC in all dimensions.
    &END CELL

    &COORD
    UNIT angstrom
H -0.567712 -0.469646 -0.645913
H 0.626116 -0.687796 0.308193
O 0 0 0
(...)
K 2.1035 2.1035 4.2735
Cl 4.1035 4.1035 2.6035                 ###
H 3.73881 3.10388 5.46104
H 3.65742 2.89924 6.989
O 3.1035 3.1035 6.207
(...)
    &END COORD

    &KIND O
      BASIS_SET DZVP-MOLOPT-GTH-q6
      POTENTIAL GTH-PBE-q6
    &END KIND
    &KIND H
      BASIS_SET DZVP-MOLOPT-GTH-q1
      POTENTIAL GTH-PBE-q1
    &END KIND
     &KIND K
      BASIS_SET DZVP-MOLOPT-SR-GTH-q9
      POTENTIAL GTH-PBE-q9
    &END KIND
    &KIND Cl
      BASIS_SET DZVP-MOLOPT-GTH-q7
      POTENTIAL GTH-PBE-q7
    &END KIND





  &END SUBSYS

 &DFT

    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  GTH_POTENTIALS
    
!    SPIN_POLARIZED ! Do spin-polarized calculation

    &POISSON
       PERIODIC XYZ
    &END POISSON

    &QS
      METHOD GPW
      EPS_DEFAULT 1.0E-10   ! Set various epsilons for QS to values that will lead
           ! to energy correct up to 1e-10.
    &END QS

    &MGRID
      CUTOFF 400    ! This is Ecut of eq. 39 in VandeVondele (2005), i.e., plane-wave cutoff
             ! that determines size of finest grid (see caption of Fig. 1). Cutoffs for
   ! the subsequent, coarser grid levels are given by eq. 39.
      NGRIDS 4      ! This is N of eq. 39 in VandeVondele (2005), i.e., number of grids used.
      REL_CUTOFF 40 ! This controls the grid level onto which Gaussians will be mapped.
    &END MGRID

    &XC

      &XC_FUNCTIONAL
&PBE
          PARAMETRIZATION ORIG
&END PBE
      &END XC_FUNCTIONAL

      &VDW_POTENTIAL

         POTENTIAL_TYPE PAIR_POTENTIAL

         &PAIR_POTENTIAL
            TYPE DFTD3
            REFERENCE_FUNCTIONAL PBE
            CALCULATE_C9_TERM .TRUE.
            PARAMETER_FILE_NAME dftd3.dat
            R_CUTOFF 15.0
         &END PAIR_POTENTIAL

      &END VDW_POTENTIAL

    &END XC

    &SCF

      SCF_GUESS RESTART ! Use data from previous run as initial guess for wavefunction.
      EPS_SCF 1.0E-6 ! Threshold for converged total energy.
      MAX_SCF 300 ! Maximum number of SCF iterations performed.

      &OT
        PRECONDITIONER NONE ! This should be stable with respect to the "Cholesky errors"
      &END OT

       &PRINT
&RESTART ON

BACKUP_COPIES 1

      &EACH
MD 1
               &END EACH

ADD_LAST NUMERIC

&END RESTART

       &END PRINT

    &END SCF

    &PRINT
      &E_DENSITY_CUBE

        STRIDE 1 1 1

        &EACH    
MD 99999999
        &END EACH

        ADD_LAST NUMERIC
     &END E_DENSITY_CUBE

      &PDOS
            COMPONENTS .FALSE.
            NLUMO = -1
            FILENAME dosfile
            LOG_PRINT_KEY TRUE

   &EACH
MD 99999999
   &END EACH

   ADD_LAST NUMERIC
      &END PDOS

    &END PRINT

 &END DFT

&END FORCE_EVAL

&MOTION

&MD
ENSEMBLE NVT
STEPS 10000
TEMPERATURE 300.0 ! K
TIMESTEP 0.5 ! fs
&THERMOSTAT

REGION GLOBAL
TYPE NOSE

&NOSE
LENGTH 3 ! Length of Nose-Hoover chain
TIMECON 20.0 ! Period of typical vibrational motion in system in fs
&END NOSE
&END THERMOSTAT
&END MD

        &PRINT

&RESTART

&EACH
     MD 1
&END EACH

ADD_LAST NUMERIC

&END RESTART

                &TRAJECTORY ON
                           ADD_LAST NUMERIC
                           FILENAME trajectory
                &END TRAJECTORY

        &END PRINT

&END MOTION


hut...@chem.uzh.ch

unread,
Jul 1, 2016, 8:10:52 AM7/1/16
to cp...@googlegroups.com
Hi

no this shouldn't be, but without more information I will have to guess.
You could also have a look at the timings at the end of the output to
see if some routines got slower or if all parts of the run were affected.

Two things to consider:

1) Use REFERENCE_C9_TERM TRUE in order to reduce the time for vdW in MD.
2) 30 SCF iterations in MD for such a simple system is pointing to a problem
with your setup.

regards

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----To: cp2k <cp...@googlegroups.com>
From: Simiam Ghan
Sent by: cp...@googlegroups.com
Date: 07/01/2016 01:38PM
Subject: [CP2K:7883] SCF, MD run-time verses atomic species.

Dear all,I am running NVT MD with Quickstep on a box of water with 64 molecules.   If I replace a water molecule with a KCl ion pair, i observe that the MD and SCF step times more than double on my setup.  (MD from ~20 sec to ~46 sec). SCF iterations are converging (except the very first MD step) in around 30 steps in both cases but each SCF now takes over twice as long as before.   Is there an explanation of why such a 'small' change in system could double the run time?  Is KCl really so heavy to calculate compared to H2O?  The number of electrons in the system increases by 8, from 512 to 520.  Below my KCl input file. 
Greetings,Simiam


&GLOBAL
  PROJECT H2O_KCl  RUN_TYPE MD  PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
  METHOD Quickstep ! GPW method.
  &SUBSYS                                       ! A subsystem: coordinates, topology, molecules and cell.
    &CELL                                       ! Supercell setup.      ABC [angstrom] 12.414 12.414 12.414 ! Using 64 H2O molecules, we thus get a density of 1g/cm^3.      PERIODIC XYZ           ! Use PBC in all dimensions.    &END CELL
    &COORD    UNIT angstromH -0.567712 -0.469646 -0.645913H 0.626116 -0.687796 0.308193O 0 0 0(...)K 2.1035 2.1035 4.2735Cl 4.1035 4.1035 2.6035                 ###H 3.73881 3.10388 5.46104
H 3.65742 2.89924 6.989O 3.1035 3.1035 6.207(...)    &END COORD
--

You received this message because you are subscribed to the Google Groups "cp2k" group.

To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.

To post to this group, send email to cp...@googlegroups.com.

Visit this group at https://groups.google.com/group/cp2k.

For more options, visit https://groups.google.com/d/optout.


Samuel Andermatt

unread,
Jul 4, 2016, 5:08:52 AM7/4/16
to cp2k
Could you provide the timing sections at the end of your calculations?

Simiam Ghan

unread,
Aug 17, 2016, 7:08:12 PM8/17/16
to cp2k
Hello Juerg,
Thank you for your reply.  I added REFERENCE_C9_TERM .TRUE.  to the vdw section as you suggested and found that my MD simulation became twice (!) as fast.   That is, time per MD step went from ~60 sec to ~30 sec for box of water with 64 molecules and K,Cl ions.  I made a quick comparison of Energy_Force results with and without this flag and indeed forces and energies do not change significantly.   This is great news.  Why is this not a default?  What's the catch?   Also, would you say this is now a correct way to declare a PBE-D3 setup?


&VDW_POTENTIAL
         POTENTIAL_TYPE  PAIR_POTENTIAL
         &PAIR_POTENTIAL
           R_CUTOFF     1.5000000000000005E+01
           TYPE  DFTD3
           PARAMETER_FILE_NAME dftd3.dat
           REFERENCE_FUNCTIONAL PBE
           CALCULATE_C9_TERM  T
           REFERENCE_C9_TERM  T
         &END PAIR_POTENTIAL
       &END VDW_POTENTIAL

Finally, I mentioned earlier that 30 scf steps were needed to converge the same box of water with ions.   I was then using EPS_SCF = 1E-06.   Reducing this to EPS_SCF = 1E-05 brought the number down to 10-20 scf necessary.  Does that still sound high?   Is the definition of EPS_SCF documented somewhere, as it is apparently not just Energy convergence. 

Cheers,
Simiam Ghan

hut...@chem.uzh.ch

unread,
Aug 18, 2016, 8:05:57 AM8/18/16
to cp...@googlegroups.com
Hi

The PBE-D3 setup looks fine.

Yes, 10-20 SCF iterations during an MD are on the very high end for
such a well behaved system. I would try the following setup (in case you
don't use it already)

&QS
....
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 4
&END QS

&OT
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
&END OT

regards

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----To: cp2k <cp...@googlegroups.com>
From: Simiam Ghan
Sent by: cp...@googlegroups.com
Date: 08/18/2016 01:08AM
Subject: Re: [CP2K:8075] SCF, MD run-time verses atomic species.

Hello Juerg,Thank you for your reply.  I added REFERENCE_C9_TERM .TRUE.  to the vdw section as you suggested and found that my MD simulation became twice (!) as fast.   That is, time per MD step went from ~60 sec to ~30 sec for box of water with 64 molecules and K,Cl ions.  I made a quick comparison of Energy_Force results with and without this flag and indeed forces and energies do not change significantly.   This is great news.  Why is this not a default?  What's the catch?   Also, would you say this is now a correct way to declare a PBE-D3 setup?

&VDW_POTENTIAL         POTENTIAL_TYPE  PAIR_POTENTIAL         &PAIR_POTENTIAL           R_CUTOFF     1.5000000000000005E+01           TYPE  DFTD3           PARAMETER_FILE_NAME dftd3.dat           REFERENCE_FUNCTIONAL PBE           CALCULATE_C9_TERM  T           REFERENCE_C9_TERM  T         &END PAIR_POTENTIAL       &END VDW_POTENTIAL
Finally, I mentioned earlier that 30 scf steps were needed to converge the same box of water with ions.   I was then using EPS_SCF = 1E-06.   Reducing this to EPS_SCF = 1E-05 brought the number down to 10-20 scf necessary.  Does that still sound high?   Is the definition of EPS_SCF documented somewhere, as it is apparently not just Energy convergence. 
Cheers,Simiam Ghan

Simiam Ghan

unread,
Aug 18, 2016, 1:30:55 PM8/18/16
to cp2k
Hi,
Thank you for the tips!  I am not using EXTRAPOLATION or PRECONDITIONER currently.   If I include them now will they affect the accuracy significantly?  I am in the middle of long MD trajectories of TiO2 crystal with water and KCl ions.  Also running the water box with ions.   My current setup is below (same for both systems).    The TiO2 system uses 1-80 scf steps for each MD step. Would you recommend these changes for that system also?  Best regards.    


 &GLOBAL
   PRINT_LEVEL  MEDIUM
   PROJECT_NAME brookite
   RUN_TYPE  MD
 &END GLOBAL
 &MOTION
   &MD
     ENSEMBLE  NVT
     STEPS  250
     TIMESTEP     4.9999999999999989E-01
     STEP_START_VAL  3751
     TIME_START_VAL     4.3234999999995443E+03
     ECONS_START_VAL    -1.5206763646667368E+04
     TEMPERATURE     3.0000000000000000E+02
     &THERMOSTAT
       TYPE  NOSE
       REGION  GLOBAL
       &NOSE
         LENGTH  3
         TIMECON     1.9999999999999993E+01
         &COORD
               4.2157834584794962E-01   -3.8977051434234506E+00    3.0959934319656378E+02
         &END COORD
         &VELOCITY
               5.7927878733942134E-05    7.8527005854180164E-04    5.0048257140974906E-04
         &END VELOCITY
         &MASS
               1.4749962622867215E+06    6.4949196930282767E+02    6.4949196930282767E+02
         &END MASS
         &FORCE
               4.2152207720701857E-08    6.1579061284564603E-06   -8.4610128141038230E-07
         &END FORCE
       &END NOSE
     &END THERMOSTAT
     &AVERAGES  T
       &RESTART_AVERAGES
         ITIMES_START  1
         AVECPU     9.1875550803485453E+01
         AVEHUGONIOT     0.0000000000000000E+00
         AVETEMP_BARO     0.0000000000000000E+00
         AVEPOT    -1.5208725197223592E+04
         AVEKIN     1.0826866556597368E+00
         AVETEMP     3.0108762437707566E+02
         AVEKIN_QM     0.0000000000000000E+00
         AVETEMP_QM     0.0000000000000000E+00
         AVEVOL     5.5116376409126919E+04
         AVECELL_A     2.7174261790895553E+01
         AVECELL_B     6.9659084710430690E+01
         AVECELL_C     2.9116900255502028E+01
         AVEALPHA     9.0000000000000071E+01
         AVEBETA     9.0000000000000071E+01
         AVEGAMMA     9.0000000000000071E+01
         AVE_ECONS     1.6039009887889247E+02
         AVE_PRESS     0.0000000000000000E+00
         AVE_PXX     0.0000000000000000E+00
       &END RESTART_AVERAGES
     &END AVERAGES
   &END MD
   &CONSTRAINT
     &FIXED_ATOMS
       COMPONENTS_TO_FIX  XYZ
       LIST  148 172 196 150 174 198 64 88 112 \
        259 270 281 221 234 247 149 173 197 289 \
        313 337 260 271 282 225 238 251 155 179 \
        203 163 187 211 164 188 212 290 314 338 \
        71 95 119 266 277 288 265 276 287
     &END FIXED_ATOMS
   &END CONSTRAINT
   &PRINT
     &TRAJECTORY  ON
       ADD_LAST  NUMERIC
       FILENAME trajectory
     &END TRAJECTORY
     &RESTART  SILENT
       ADD_LAST  NUMERIC
       &EACH
         MD  1
       &END EACH
     &END RESTART
   &END PRINT
 &END MOTION
 &FORCE_EVAL
   METHOD  QS
   &DFT
     BASIS_SET_FILE_NAME BASIS_MOLOPT
     POTENTIAL_FILE_NAME GTH_POTENTIALS
     &SCF
       MAX_SCF  300
       EPS_SCF     1.0000000000000001E-05
       SCF_GUESS  RESTART
       &OT  T
         PRECONDITIONER  NONE
       &END OT
       &PRINT
         &RESTART  ON
           ADD_LAST  NUMERIC
           BACKUP_COPIES  1
           &EACH
             MD  1
           &END EACH
         &END RESTART
       &END PRINT
     &END SCF
     &QS
       EPS_DEFAULT     1.0000000000000000E-10
       METHOD  GPW
     &END QS
     &MGRID
       NGRIDS  4
       CUTOFF     4.0000000000000000E+02
       REL_CUTOFF     4.0000000000000000E+01
     &END MGRID
     &XC
       DENSITY_CUTOFF     1.0000000000000000E-10
       GRADIENT_CUTOFF     1.0000000000000000E-10
       TAU_CUTOFF     1.0000000000000000E-10
       &XC_FUNCTIONAL  NO_SHORTCUT
         &PBE  T
           PARAMETRIZATION  ORIG
         &END PBE
       &END XC_FUNCTIONAL
       &VDW_POTENTIAL
         POTENTIAL_TYPE  PAIR_POTENTIAL
         &PAIR_POTENTIAL
           R_CUTOFF     1.5000000000000005E+01
           TYPE  DFTD3
           PARAMETER_FILE_NAME dftd3.dat
           REFERENCE_FUNCTIONAL PBE
           CALCULATE_C9_TERM  T
           REFERENCE_C9_TERM  T
         &END PAIR_POTENTIAL
       &END VDW_POTENTIAL
     &END XC
     &POISSON
       PERIODIC  XYZ
     &END POISSON
     &PRINT
       &E_DENSITY_CUBE  SILENT
         ADD_LAST  NUMERIC
         STRIDE  1 1 1
         &EACH
           MD  99999999
         &END EACH
       &END E_DENSITY_CUBE
       &PDOS  SILENT
         ADD_LAST  NUMERIC
         FILENAME dosfile
         LOG_PRINT_KEY  T
         COMPONENTS  F
         NLUMO  -1
         &EACH
           MD  99999999
         &END EACH

hut...@chem.uzh.ch

unread,
Aug 19, 2016, 3:11:48 AM8/19/16
to cp...@googlegroups.com
Hi

Setting an apropriate preconditioner and the DIIS minimizer will
speed up convergence and will not affect accuracy.
The ASPC order=4 setting gives you better initial wavefunctions
(=faster convergence) and better energy conservation (=stability of MD).
Again accuracy will be the same or better.

regards

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----To: cp2k <cp...@googlegroups.com>
From: Simiam Ghan
Sent by: cp...@googlegroups.com
Date: 08/18/2016 07:31PM
Subject: Re: [CP2K:8079] SCF, MD run-time verses atomic species.

Hi,Thank you for the tips!  I am not using EXTRAPOLATION or PRECONDITIONER currently.   If I include them now will they affect the accuracy significantly?  I am in the middle of long MD trajectories of TiO2 crystal with water and KCl ions.  Also running the water box with ions.   My current setup is below (same for both systems).    The TiO2 system uses 1-80 scf steps for each MD step. Would you recommend these changes for that system also?  Best regards.    

Simiam Ghan

unread,
Aug 21, 2016, 1:36:13 PM8/21/16
to cp2k
Greetings Juerg, 
I attached a picture showing the effect of REFERENCE_C9_TERM and EXTRAPOLATION and PRECONDITIONER on my MD simulation.  Together they have reduced time per MD step from 75-200sec to 36sec and SCF count from 1-80 to 6.    Now, as this seems too good to be true I did an ENERGY_FORCE calculation on a snapshot of the system with and without REFERENCE_C9_TERM: the difference in total energy is  0.02meV/atom.  With/without the EXTRAPOLATION and PRECONDITIONER which you gave me the total energy difference is 1meV/atom.  This difference is reduced if i tighten EPS_SCF from 1E-5 to 1E-6. I have not checked if these small changes translate to better/worse accuracy. 

Could you elaborate on what REFERENCE_C9_TERM does: how does it double speed while maintaining accuracy? Why is it not a default? 
Also , could you explain the definition of EPS_SCF, for example what are its units?

Many thanks,
Simiam

P.S.  EXTRAPOLATION and PRECONDITIONER also put a stop to drifting in the Cons.Quantity of the .ener file, thank you. 
SCf_performance.pdf

hut...@chem.uzh.ch

unread,
Aug 23, 2016, 3:43:14 AM8/23/16
to cp...@googlegroups.com
Hi

REFERENCE_C9_TERM True

keeps the C6 terms used for the calculation of the C9 term constant.
This considerably simplifies the calculation of forces.

EPS_SCF in OT is MAXVAL(dE/dCai - lambda Sab*Cbi) [the 'gradient'] in
atomic units.

regards

Juerg
--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie C                FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----To: cp2k <cp...@googlegroups.com>
From: Simiam Ghan
Sent by: cp...@googlegroups.com
Date: 08/21/2016 07:36PM
Subject: Re: [CP2K:8085] SCF, MD run-time verses atomic species.
[attachment "SCf_performance.pdf" removed by Jürg Hutter/at/UZH]
Reply all
Reply to author
Forward
0 new messages