Force calculation in CELL_OPT vs MD

[martin....@gmail.com]

Hello,

This may well just be a gap in my knowledge regarding how forces are calculated in CP2K for geometry optimisations vs molecular dynamics. In such a case, I would still appreciate any insight as well.

In short, for the same system with virtually the same SCF and other settings, force calculations during MD take over 4 times as long as the same step during CELL_OPT. I don't understand why this is happening and would appreciate any insight or solution.

Details: I have started by optimising the geometry of this system (relatively tight convergence criteria, no issues) prior to running any MD. My issue comes during the MD, where each step takes significantly longer (average of ~30 seconds per opt. step vs ~90 seconds per md step) than i would expect based on the corresponding CELL_OPT timings. Most of this seems to be taken up by force calculations (timings in the output file suggest ~20 seconds for the SCF process and finding the electron density and ~60-70 seconds for the force calculation, position update etc.). Previously, similar MD calcs on a cluster model of this system (but containing more atoms in a larger unit cell) required 35-45 s per time step. This is despite using the same numbers of cores, same basis sets, same cutoffs, and same settings where relevant.

Background: I am working on a certain porous material (unit cell 15x15x20) with regards to catalytic processes which can take place inside the pores. CP2K 6.1 or 7.1 (issue persists with both), 600/50 Ry cutoff and rel cutoff, PBE-d3 functional with DZVP-MOLOPT basis sets, BFGS optimiser for CELL_OPT,  GTH pseudopotentials. CSVR thermostat with 0.5 fs time step and 1 fs time constant.

Best regards,

Martin

[Krack Matthias (PSI)]

Hi Martin

 

The probability of getting any feedback on this forum would increase, if you would provide the input files and the timings in the output files of the CELL_OPT and MD runs.

 

Best

 

Matthias

--
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 view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/3964b36a-2b55-4013-9573-e05ad5f58268n%40googlegroups.com.

[martin....@gmail.com]

Hello,

Below are the inputs and timings for the cell opt and 20 subsequent md steps:

CELL_OPT:

&GLOBAL
    PROJECT uio66-half
    RUN_TYPE CELL_OPT
    PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
    METHOD QS
    &SUBSYS
        &CELL
            ABC 14.8 14.8 20.910
            ALPHA_BETA_GAMMA 90 90 90
        &END CELL
        &TOPOLOGY
            COORD_FILE_FORMAT CIF
            COORD_FILE_NAME uio66-pristine-half-cell.cif
        &END TOPOLOGY
        &KIND O
            BASIS_SET DZVP-MOLOPT-GTH
            POTENTIAL GTH-PBE-q6
        &END KIND
        &KIND Zr
            BASIS_SET DZVP-MOLOPT-SR-GTH
            POTENTIAL GTH-PBE-q12
        &END KIND
        &KIND C
            BASIS_SET DZVP-MOLOPT-GTH
            POTENTIAL GTH-PBE-q4
        &END KIND
        &KIND H
            BASIS_SET DZVP-MOLOPT-GTH
            POTENTIAL GTH-PBE-q1
        &END KIND
    &END SUBSYS
    &DFT
        BASIS_SET_FILE_NAME BASIS_MOLOPT
        POTENTIAL_FILE_NAME GTH_POTENTIALS
        &QS
            EPS_DEFAULT 1.0E-12
        &END QS
        &MGRID
            CUTOFF 600
            NGRIDS 5
            REL_CUTOFF 50
        &END MGRID
        &SCF
            SCF_GUESS ATOMIC
            EPS_SCF 1.0E-08
            MAX_SCF 30
            &OT
                MINIMIZER DIIS
                                PRECONDITIONER FULL_KINETIC
                LINESEARCH 2PNT
                        &END OT
            &OUTER_SCF
                EPS_SCF 1.0E-08
                MAX_SCF 25
            &END OUTER_SCF
        &END SCF
        &XC
            &XC_FUNCTIONAL PBE
            &END XC_FUNCTIONAL
            &vdW_POTENTIAL
                DISPERSION_FUNCTIONAL PAIR_POTENTIAL
                &PAIR_POTENTIAL
                    TYPE DFTD3
                    CALCULATE_C9_TERM TRUE
                    PARAMETER_FILE_NAME dftd3.dat
                    REFERENCE_FUNCTIONAL PBE
                &END PAIR_POTENTIAL
            &END vdW_POTENTIAL
        &END XC
    &END DFT
    STRESS_TENSOR ANALYTICAL
&END FORCE_EVAL
&MOTION
    &CELL_OPT
                KEEP_ANGLES TRUE
        TYPE DIRECT_CELL_OPT
                MAX_DR 4.5E-04
                MAX_FORCE 1.0E-04
                RMS_DR 4.5E-04
                RMS_FORCE 4.5E-04
        MAX_ITER 400
        OPTIMIZER BFGS
    &END CELL_OPT
&END MOTION

MD

&GLOBAL
   PRINT_LEVEL  LOW
   PROJECT_NAME uio66-defect-md
   RUN_TYPE  MD
 &END GLOBAL
 &MOTION
   &MD
     ENSEMBLE  NVT
     STEPS  20
     TIMESTEP     4.9999999999999989E-01
     STEP_START_VAL  20
     TIME_START_VAL     1.0000000000000002E+01
     ECONS_START_VAL    -2.2805855491897287E+03
     TEMPERATURE     3.0000000000000000E+02
     &THERMOSTAT
       TYPE  CSVR
       &CSVR
         TIMECON     9.9999999999999978E-01
         &THERMOSTAT_ENERGY
              -3.1503507886676752E-01
         &END THERMOSTAT_ENERGY
         &RNG_INIT
Wiener process for Thermostat # 1        1 F T F  -5.2769235486593846E-01        1600409560.0         760341491.0        3744745813.0        3560625215.0         548769459.0        3217800586.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0             12345.0
         &END RNG_INIT
       &END CSVR
     &END THERMOSTAT
     &AVERAGES  T
       &RESTART_AVERAGES
         ITIMES_START  1
         AVECPU     9.3657589755114145E+01
         AVEHUGONIOT     0.0000000000000000E+00
         AVETEMP_BARO     0.0000000000000000E+00
         AVEPOT    -2.1740421539898202E+03
         AVEKIN     3.2332755503992078E-01
         AVETEMP     2.9984917723950269E+02
         AVEKIN_QM     0.0000000000000000E+00
         AVETEMP_QM     0.0000000000000000E+00
         AVEVOL     3.0842651782127203E+04
         AVECELL_A     2.7971726218973288E+01
         AVECELL_B     2.7971726218973288E+01
         AVECELL_C     3.9419687131994479E+01
         AVEALPHA     9.0000000000000000E+01
         AVEBETA     9.0000000000000000E+01
         AVEGAMMA     9.0000000000000000E+01
         AVE_ECONS     1.4016569135616819E+05
         AVE_PRESS    -1.6562449007126843E+03
         AVE_PXX    -3.4950333563135082E+03
       &END RESTART_AVERAGES
     &END AVERAGES
   &END MD
   &PRINT
     &TRAJECTORY  SILENT
       &EACH
         MD  1
       &END EACH
     &END TRAJECTORY
     &VELOCITIES  SILENT
       &EACH
         MD  500
       &END EACH
     &END VELOCITIES
     &RESTART  SILENT
       &EACH
         MD  1
       &END EACH
     &END RESTART
     &RESTART_HISTORY  SILENT
       &EACH
         MD  50
       &END EACH
     &END RESTART_HISTORY
   &END PRINT
 &END MOTION
 &FORCE_EVAL
   METHOD  QS
   STRESS_TENSOR  ANALYTICAL
   &DFT
     BASIS_SET_FILE_NAME BASIS_MOLOPT
     POTENTIAL_FILE_NAME POTENTIAL
     &SCF
       MAX_SCF  30
       EPS_SCF     1.0000000000000000E-08
       SCF_GUESS  RESTART
       &OT  T
         MINIMIZER  DIIS
         LINESEARCH  2PNT
         PRECONDITIONER  FULL_ALL
       &END OT
       &OUTER_SCF  T
         EPS_SCF     1.0000000000000000E-08
         MAX_SCF  100
       &END OUTER_SCF
     &END SCF
     &QS
       EPS_DEFAULT     1.0000000000000002E-12
       EPS_PGF_ORB     1.0000000000000001E-18
     &END QS
     &MGRID
       NGRIDS  5
       CUTOFF     6.0000000000000000E+02
       REL_CUTOFF     5.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
         &END PBE
       &END XC_FUNCTIONAL
       &VDW_POTENTIAL
         POTENTIAL_TYPE  PAIR_POTENTIAL
         &PAIR_POTENTIAL
           TYPE  DFTD3
           PARAMETER_FILE_NAME dftd3.dat
           REFERENCE_FUNCTIONAL PBE
           CALCULATE_C9_TERM  T
         &END PAIR_POTENTIAL
       &END VDW_POTENTIAL
     &END XC
   &END DFT
   &SUBSYS
     &CELL
       A     1.4802000000000001E+01    0.0000000000000000E+00    0.0000000000000000E+00
       B     0.0000000000000000E+00    1.4802000000000001E+01    0.0000000000000000E+00
       C     0.0000000000000000E+00    0.0000000000000000E+00    2.0859999999999999E+01
       MULTIPLE_UNIT_CELL  1 1 1
     &END CELL
     &KIND H
       BASIS_SET DZVP-MOLOPT-GTH
       POTENTIAL GTH-PBE-q1
     &END KIND
     &KIND O
       BASIS_SET DZVP-MOLOPT-GTH
       POTENTIAL GTH-PBE-q6
     &END KIND
     &KIND C
       BASIS_SET DZVP-MOLOPT-GTH
       POTENTIAL GTH-PBE-q4
     &END KIND
     &KIND Zr
       BASIS_SET DZVP-MOLOPT-SR-GTH
       POTENTIAL GTH-PBE-q12
     &END KIND
     &TOPOLOGY

       COORD_FILE_FORMAT XYZ

       COORD_FILE_NAME half-cell-optimised.xyz

       NUMBER_OF_ATOMS  228
       MULTIPLE_UNIT_CELL  1 1 1
     &END TOPOLOGY
   &END SUBSYS
 &END FORCE_EVAL

The latter is a restart file from which ive removed coordinates and velocities to reduce cluttering. The preconditioner changes from full_kinetic to full_all for md, but using full_kinetic for both leads to no major difference. The timings are attached as text files to preserve sensible formatting.

Best regards,

Martin

[Krack Matthias (PSI)]

Hi Martin

 

The timings for the core Hamiltonian forces, especially the nonlocal PP part (pnnl), are quite large in the MD run. This is most likely due to the tight EPS_PGF_ORB parameter (1.0E-18). Remove that parameter in the MD input and just rely on EPS_DEFAULT like in the CEL_OPT run. The *_CUTOFF parameters in the &XC section are also dispensable.

 

HTH

 

Matthias

--

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 view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e9f76caa-dbce-49a3-a778-17912e6c2915n%40googlegroups.com.

[martin....@gmail.com]

Hello Matthias,

Your suggestion seems to have done the trick. This hasn't caused me any problems for the other system i described, but thinking about it now, the distribution of atoms there is such that the tight eps convergence for overlap elements is not so much of an issue. Thank you for your help.

Best regards,

Martin