SCCS Solvation Model Implementation - Molecule Radii Explodes

216 views
Skip to first unread message
Assigned to dr.fen...@gmail.com by me

Dev Rana

unread,
May 11, 2020, 11:24:38 AM5/11/20
to cp2k
Hello,

I'm trying to simulate a single carbon atom within an aluminum block using a solvation model for the aluminum. I can do this in a normal MD simulation, however to reduce computation time/cost I'd like to use a solvation model to reduce system size. Using the QS regtest-3 files as a starting point, I've created my own input file for my, now 5 atom (Al4C) system. I am finding that the bond lengths between my carbon and aluminum atoms go from something reasonable to over 5 A, and that GEO_OPT is still not complete. Any suggestions would be appreciated. I've attached the input, output, and trajectory files for your review. I've also paste the input below.

Best Regards,
Dev

&GLOBAL
 PROJECT SolvationAl4C
 RUN_TYPE geo_opt
 PRINT_LEVEL medium
&END GLOBAL


&FORCE_EVAL
 METHOD Quickstep
 # STRESS_TENSOR analytical
 &DFT
    BASIS_SET_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/BASIS_MOLOPT
    POTENTIAL_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/POTENTIAL
  
  &QS
    METHOD GPW
    EPS_DEFAULT 1.0E-10
    EXTRAPOLATION ASPC
  &END QS
  
  &MGRID
    NGRIDS 4
    CUTOFF [Ry] 600
    REL_CUTOFF [Ry] 250
   &RS_GRID
    DISTRIBUTION_TYPE distributed
   &END RS_GRID
  &END MGRID

  &SCCS ON
   ALPHA [N*m^-1] 0.0
   BETA [GPa] 1.3
   DELTA_RHO 2.0E-5
   DERIVATIVE_METHOD fft
   DIELECTRIC_CONSTANT 80
   EPS_SCCS 1.0E-3
   GAMMA [mN/m] 0.0
   EPS_SCF 0.0
   MAX_ITER 100
!   METHOD Andreussi
   METHOD Fattebert-Gygi
   MIXING 0.6
   &ANDREUSSI
    RHO_MAX 0.001
    RHO_MIN 0.0001
   &END ANDREUSSI
   &FATTEBERT-GYGI
    BETA 1.3
    RHO_ZERO 0.00078
   &END FATTEBERT-GYGI
  &END SCCS
  
  &SCF
    SCF_GUESS ATOMIC
    MAX_SCF 300
    EPS_SCF 1.0E-3
!    ADDED_MOS 200
!    &DIAGONALIZATION
!      ALGORITHM STANDARD
!    &END DIAGONALIZATION
!    &SMEAR TRUE
!     METHOD FERMI_DIRAC
!     ELECTRONIC_TEMPERATURE 300
!    &END SMEAR
!    &MIXING
!      METHOD BROYDEN_MIXING
!      ALPHA 0.1
!      BETA 1.0
!      NBROYDEN 8
!    &END MIXING
    &OT on
     MINIMIZER DIIS
     PRECONDITIONER FULL_SINGLE_INVERSE
    &END OT
  &END SCF
  
  &XC
   &XC_FUNCTIONAL PBE
   &END XC_FUNCTIONAL
  &END XC
  
  &PRINT
    &PDOS
      &EACH
        GEO_OPT 1
      &END EACH
      APPEND TRUE
    &END PDOS
    &SCCS ON
      &EACH
       GEO_OPT 10
      &END EACH
     &DENSITY_GRADIENT on
      &EACH
       GEO_OPT 10
      &END EACH
      APPEND
     &END DENSITY_GRADIENT
     &DIELECTRIC_FUNCTION on
      &EACH
       GEO_OPT 10
      &END EACH
      APPEND
     &END DIELECTRIC_FUNCTION
     &POLARISATION_POTENTIAL on
      &EACH
       GEO_OPT 10
      &END EACH
      APPEND
     &END POLARISATION_POTENTIAL
    &END SCCS
  &END PRINT
  
 &END DFT
 
 &SUBSYS
    &CELL
      ABC [angstrom] 10 10 10
      PERIODIC XYZ
      MULTIPLE_UNIT_CELL 1 1 1
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME Al4C.xyz
      COORD_FILE_FORMAT XYZ
      MULTIPLE_UNIT_CELL 1 1 1
    &END
    &KIND C
      ELEMENT C
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q4
    &END KIND
    &KIND Al
      ELEMENT Al
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q3
    &END KIND
 &END SUBSYS

&END FORCE_EVAL

&MOTION
 &GEO_OPT
  OPTIMIZER BFGS
  MAX_ITER  100
 &END GEO_OPT
 &PRINT
  &STRUCTURE_DATA
    DISTANCE 1 2
    DISTANCE 1 3
    DISTANCE 1 4
    DISTANCE 1 5
    DISTANCE 2 3
    DISTANCE 2 4
    DISTANCE 2 5
    DISTANCE 3 4
    DISTANCE 3 5
    DISTANCE 4 5
    ANGLE 1 2 3
    ANGLE 1 2 4
    ANGLE 1 2 5
    ANGLE 1 3 2
    ANGLE 1 3 4
    ANGLE 1 3 5
    ANGLE 1 4 2
    ANGLE 1 4 3
    ANGLE 1 4 5
    ANGLE 1 5 2
    ANGLE 1 5 3
    ANGLE 1 5 4
    &EACH
      GEO_OPT 1
    &END EACH
  &END STRUCTURE_DATA
  &TRAJECTORY
    &EACH
      GEO_OPT 1
    &END EACH
  &END TRAJECTORY
  &VELOCITIES
    &EACH
      GEO_OPT 1
    &END EACH
  &END VELOCITIES
  &FORCES
    &EACH
      GEO_OPT 1
    &END EACH
  &END FORCES
  &CELL
    &EACH
      GEO_OPT 1
    &END EACH
  &END CELL
  &RESTART
    &EACH
      GEO_OPT 1
    &END EACH
  &END RESTART
 &END PRINT
&END MOTION

job.in
output.out
SolvationAl4C-pos-1.xyz

Krack Matthias (PSI)

unread,
May 11, 2020, 4:56:51 PM5/11/20
to cp...@googlegroups.com

Dear Dev

 

an EPS_SCF threshold value of 0.001 for the SCF convergence is by far too large which results in unconverged energies and forces screwing up your system. Likewise, I suggest a smaller EPS_DEFAULT value of at least 1.0E-12 or less. It is also highly recommended to use an OUTER_SCF with OT.

 

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/b2a95387-2350-4f1c-8eb9-32976ddffbd0%40googlegroups.com.

Dev Rana

unread,
May 13, 2020, 3:58:32 PM5/13/20
to cp2k
Thanks Matthias!

I reran using your suggestions and it looks like SCF convergence does not complete. I also tried to turn off OT and use Diagonalization with Outer SCF too. But it all results in similar nonconvergence. I ran a GEO_OPT without SCCS and it does indeed converge. However upon adding SCCS, it just computes forever with each SCF step taking over 300 seconds each. I've attached my files for your reference and my input below.

What do you suggest?

&GLOBAL
 PROJECT SCCS
 RUN_TYPE geo_opt
 PRINT_LEVEL low
 SAVE_MEM TRUE
&END GLOBAL


&FORCE_EVAL
 METHOD Quickstep
 # STRESS_TENSOR analytical
 &DFT
    BASIS_SET_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/BASIS_MOLOPT
    POTENTIAL_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/POTENTIAL
  
  &QS
    METHOD GPW
    EPS_DEFAULT 1.0E-12
    EXTRAPOLATION ASPC
  &END QS
  
  &MGRID
    NGRIDS 5
    CUTOFF [Ry] 600
    REL_CUTOFF [Ry] 250
  &END MGRID

  &SCCS
   ALPHA [N*m^-1] 0.0
   BETA [GPa] 0.0
   DELTA_RHO 2.0E-5
   DERIVATIVE_METHOD fft
   DIELECTRIC_CONSTANT 80
   GAMMA [mN/m] 0.0
   EPS_SCCS 1.0E-6
   EPS_SCF 0.0
   MAX_ITER 100
!   METHOD Andreussi
   METHOD Fattebert-Gygi
   MIXING 0.6
!   &ANDREUSSI
!    RHO_MAX 0.001
!    RHO_MIN 0.0001
!   &END ANDREUSSI
   &FATTEBERT-GYGI
    BETA 1.3
    RHO_ZERO 0.00078
   &END FATTEBERT-GYGI
  &END SCCS
  
  &SCF
    SCF_GUESS ATOMIC
    MAX_SCF 300
    EPS_SCF 1.0E-7
    &DIAGONALIZATION ON
      ALGORITHM STANDARD
    &END DIAGONALIZATION
!    &OT on
!     MINIMIZER DIIS
!     PRECONDITIONER FULL_ALL
!     ENERGY_GAP 0.002
!    &END OT
    &OUTER_SCF
     EPS_SCF 1e-7
     MAX_SCF 30
    &END OUTER_SCF
      ABC [angstrom] 30 30 30
!      PERIODIC XYZ
!      MULTIPLE_UNIT_CELL 1 1 1
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME Al4C.xyz
      COORD_FILE_FORMAT XYZ
!      MULTIPLE_UNIT_CELL 1 1 1
    &END
    &KIND C
      ELEMENT C
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q4
    &END KIND
    &KIND Al
      ELEMENT Al
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q3
    &END KIND
 &END SUBSYS

&END FORCE_EVAL

&MOTION
 &GEO_OPT
  TYPE MINIMIZATION
  OPTIMIZER BFGS
  MAX_ITER  100
 &END GEO_OPT
 &CONSTRAINT
   &FIXED_ATOMS
     COMPONENTS_TO_FIX XYZ
     LIST 1
   &END FIXED_ATOMS
 &END CONSTRAINT
Al4C.xyz
job.in
output.out

Dev Rana

unread,
May 13, 2020, 4:40:49 PM5/13/20
to cp2k
After typing out this message. I was able to get convergence in the SCF step. The solution was to change the SCCS%EPS_SCF to 1E-7. I understood this function to mean that upon SCF convergence at the EPS_SCF value, is when the SCCS switch would turn on. Instead this serves as the same function as DFT%SCF%EPS_SCF and should be identical values.

TL;DR SCCS%EPS_SCF and DFT%SCF%EPS_SCF should be the same very small value. 

Matthias, in the regtest_3 gpw tests, the SCCS testing conducted does not accurately test out the requirements for SCCS as the EPS_SCF is set to 0.0 within that regtest file. 

Best Regards,
Dev

Dev Rana

unread,
May 13, 2020, 6:28:41 PM5/13/20
to cp2k
I take back what I previously said. The 0th geo opt step finishes quickly with SCF steps taking 20-25s each. The 2st geo opt step does not finish and each SCF step takes 5+ minutes. Here's my input file and output for that.  I'm also running this on GPU. 

What do you advise? I've attached the input output files and paste input below. The convergence values keep increases, rather than decreases. And Total Energy keep lowering. 

&GLOBAL
 PROJECT SCCS
 RUN_TYPE geo_opt
 PRINT_LEVEL low
 SAVE_MEM TRUE
&END GLOBAL


&FORCE_EVAL
 METHOD Quickstep
 # STRESS_TENSOR analytical
 &DFT
    BASIS_SET_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/BASIS_MOLOPT
    POTENTIAL_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/POTENTIAL
  
  &QS
    METHOD GPW
    EPS_DEFAULT 1.0E-12
    EXTRAPOLATION ASPC
  &END QS
  
  &MGRID
    NGRIDS 5
    CUTOFF [Ry] 600
  &END MGRID

  &SCCS
   ALPHA [N*m^-1] 0.0
   BETA [GPa] 0.0
   DELTA_RHO 2.0E-5
   DERIVATIVE_METHOD cd5
   DIELECTRIC_CONSTANT 80
   GAMMA [mN/m] 0.0
   EPS_SCCS 1.0E-6
   EPS_SCF 1.0E-7
     &DENSITY_GRADIENT off
      &EACH
       GEO_OPT 10
      &END EACH
      APPEND
     &END DENSITY_GRADIENT
     &DIELECTRIC_FUNCTION off
      &EACH
       GEO_OPT 10
      &END EACH
      APPEND
     &END DIELECTRIC_FUNCTION
     &POLARISATION_POTENTIAL off
Al4C.xyz
job.in
output_noconvergence.out

Krack Matthias (PSI)

unread,
May 13, 2020, 8:08:47 PM5/13/20
to cp...@googlegroups.com

Hi Dev

 

The settings in the regression tests are often physically not meaningful, since these test runs have to be completed after a few seconds.

SCCS%EPS_SCF delays the start of the inner SCCS convergence loop only for values greater zero. Otherwise that value is not considered and the inner SCCS cycle is activated immediately.

SCCS%EPS_SCCS, the convergence threshold of the inner SCCS loop for the polarization charge should be smaller (tighter) than 1.0E-6, at least 1.0E-8. As soon as the inner SCCS loop is activated, the timings for each SCF iteration step will increase significantly. Also the implicit solvent comes not for free. You can reduce the box size to save computer resources, because your system is small. Does the setup for your system CAl4 work/converge properly for a geo_opt run in gas phase, i.e. without SCCS? If that is the case then I would start with a smaller dielectric constant than 80 for the SCCS and see how the system evolves. The print level medium provides a detailed SCCS output.

 

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.

Dev Rana

unread,
May 15, 2020, 11:47:21 AM5/15/20
to cp2k
Hi Matthias,

Thank you again for all your help. I went a bit beyond your suggestion and changed EPS_SCF to 1E-12 and changed cell size to 10 x 10 x 10 A3. To answer your questions:

1.) I was able to run GEO_OPT without SCCS and the simulation completed. 
2.) I was able to converge and complete a SCCS GEO_OPT simulation with a dielectric constant of 1. 

I then took the dieletric constant back up to 80 (with no other changes), and the convergence fails - in fact the convergence Hartree's are trending greater than 1. 

I'm going to do a small ramp up in dielectric constant: 2, 5, 10, etc. To see what happens. But I should be able to use 80+ as the dielectric constant without issues.

What do you think? Am I missing something?

Again input, output, and XYZ attached. 

&GLOBAL
 PROJECT SCCS
 RUN_TYPE geo_opt
 PRINT_LEVEL medium
 SAVE_MEM TRUE
&END GLOBAL


&FORCE_EVAL
 METHOD Quickstep
 # STRESS_TENSOR analytical
 &DFT
    BASIS_SET_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/BASIS_MOLOPT
    POTENTIAL_FILE_NAME /shared/centos7/cp2k/cp2k-6.1.0/data/POTENTIAL
  
  &QS
    METHOD GPW
    EPS_DEFAULT 1.0E-12
    EXTRAPOLATION ASPC
  &END QS
  
  &MGRID
    NGRIDS 5
    CUTOFF [Ry] 600
  &END MGRID

  &SCCS
   ALPHA [N*m^-1] 0.0
   BETA [GPa] 0.0
   DELTA_RHO 2.0E-5
   DERIVATIVE_METHOD fft
   DIELECTRIC_CONSTANT 80
   GAMMA 0.0
   EPS_SCCS 1.0E-12
!   EPS_SCF 1.0E-10
   MAX_ITER 100
!   METHOD Andreussi
   METHOD Fattebert-Gygi
   MIXING 0.6
!   &ANDREUSSI
!    RHO_MAX 0.008
!    RHO_MIN 0.00015
!   &END ANDREUSSI
   &FATTEBERT-GYGI
    BETA 1.3
    RHO_ZERO 0.00078
   &END FATTEBERT-GYGI
  &END SCCS
  
  &SCF
    SCF_GUESS ATOMIC
    MAX_SCF 300
    EPS_SCF 1.0E-12
    &DIAGONALIZATION ON
      ALGORITHM STANDARD
    &END DIAGONALIZATION
!    &OT on
!     MINIMIZER DIIS
!     PRECONDITIONER FULL_ALL
!     ENERGY_GAP 0.002
!    &END OT
    &OUTER_SCF
     EPS_SCF 1e-12
      ABC [angstrom] 10 10 10

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

Al4C.xyz
job.in
output.out

Krack Matthias (PSI)

unread,
May 15, 2020, 1:16:45 PM5/15/20
to cp...@googlegroups.com

Hi Dev

 

The attached revised input converged after a few steps using SCCS with a dielectric constant of 80. I changed the method and adapted the rho_min and rho_max values which define the surface of the cavity and thus they are system dependent to some extent. If rho_max is chosen too large for instance, the surface might be too close to the molecule. It is often helpful to plot the cavity, e.g. with VMD, using a cube file dump of the dielectric function.

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/cc63d43c-380a-4575-ac41-5b981251d294%40googlegroups.com.

Al4C.inp

Dev Rana

unread,
May 17, 2020, 12:02:21 PM5/17/20
to cp2k
Hi Matthias,

Thank you so much! This worked. However, when I go beyond a dieletric constant of 300 i.e. 300, 400, 500, 1000+, with no other changes to the input file. The convergence starts to go beyond 1 again. Do you have a reference that might help me understand what input parameters I should be using to offset my changing conditions? I'm going to try and print a cube file as you suggested and visualize it in VMD, I'm assuming with the newer higher dieletric constants the surface might be too close to the molecule and I will have to change my rho min/max to adjust for it. however I'm not sure by how much and I hate to keep troubling you. 

Best Regards,
Dev

Krack Matthias (PSI)

unread,
May 18, 2020, 2:52:54 AM5/18/20
to cp...@googlegroups.com

Hi Dev

 

Never tried such large dielectric constants of 100-1000. Values much larger than 80 for a solvent seem to me unusual.

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/99722c46-cf0e-427f-a7f5-bd6eba017cd3%40googlegroups.com.

Dev Rana

unread,
Jun 3, 2020, 10:15:12 AM6/3/20
to cp2k
Hi Matthias,

Thank you for your help!

I was able to go to 200 without it crashing. My rationale for going above 100+ is that in the COSMO model the dielectric constant has a range that goes towards infinity because of metals. A dielectric constant of 80, is good for water but does not accurately represent an infinity model which COSMO does for metals. But I'm pushing forward with the capabilities that we have with CP2K, as I do not have access to the COSMO model. 

Best Regards,
Dev
Reply all
Reply to author
Forward
0 new messages