CP2K: How to Define a Molecule and Restrict Lennard-Jones Interactions Within It?

Planck

unread,

Feb 18, 2025, 9:06:02 AM2/18/25

to cp2k

Hi,
I am trying to simulate molecular dynamics for a hypothetical chain of carbon atoms (this is a minimum working example designed to illustrate a larger issue I am facing). Each molecule consists of 5 carbon atoms, and my simulation contains 4 such molecules.

```
20
#Built with Packmol
C 264.316076 114.220486 50.721755
C 264.392382 114.277533 50.950866
C 264.525927 114.377373 51.351839
C 264.659472 114.477213 51.752813
C 264.735777 114.534260 51.981923
C 327.956067 237.050267 49.303207
C 327.861562 237.033691 49.074379
C 327.696165 237.004681 48.673899
C 327.530769 236.975671 48.273418
C 327.436264 236.959095 48.044590
C 433.879103 142.141879 115.467902
C 433.814625 141.918476 115.381286
C 433.701780 141.527492 115.229696
C 433.588936 141.136507 115.078106
C 433.524458 140.913105 114.991490
C 234.869472 121.984410 116.856352
C 234.879704 121.977472 116.608530
C 234.897610 121.965330 116.174809
C 234.915516 121.953188 115.741088
C 234.925748 121.946251 115.493266
```
I have defined a Lennard-Jones potential between the carbon atoms. My input file looks like this:

```
&GLOBAL
PRINT_LEVEL HIGH ! Low verbosity
PROJECT MD_simulation
RUN_TYPE MD
&END GLOBAL

&MOTION
&MD
ENSEMBLE NVT ! Constant volume
TEMPERATURE 420
TIMESTEP 4.0 ! 4 fs
STEPS 100000 !

&THERMOSTAT
TYPE CSVR
&END THERMOSTAT

&PRINT
&ENERGY
&EACH
MD 25 ! 25*TIMESTEP
&END EACH
&END ENERGY
&PROGRAM_RUN_INFO
&EACH
MD 25 ! 25*TIMESTEP
&END EACH
&END PROGRAM_RUN_INFO
&END PRINT
&END MD
&PRINT
&TRAJECTORY
&EACH
MD 25 ! 25*TIMESTEP
&END EACH
&END TRAJECTORY
&END PRINT
&END MOTION

&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
DO_ELECTROSTATICS F ! Ignore electrostatics (potential from QC)
IGNORE_MISSING_CRITICAL_PARAMS T ! Ignore bonded parameters (potential from QC)

&SPLINE
EMAX_SPLINE 5000
EPS_SPLINE 1.00000000E-8 ! Increase spline accuracy
&END SPLINE

&NONBONDED

&LENNARD-JONES
atoms C C
EPSILON [kcalmol] 0.08
SIGMA [angstrom] 3.1507
RCUT [angstrom] 11.4
RMIN 1.0
&END LENNARD-JONES

&END NONBONDED

&END FORCEFIELD

&POISSON
PERIODIC NONE
&EWALD
EWALD_TYPE NONE
&END EWALD
&END POISSON

&END MM

&SUBSYS

&CELL
PERIODIC XYZ
ABC 800 800 800
&END CELL

&TOPOLOGY
COORD_FILE_FORMAT XYZ
COORD_FILE_NAME C5_20.xyz
&END TOPOLOGY

&END SUBSYS

STRESS_TENSOR NUMERICAL
&END FORCE_EVAL

```

but when I run the simulation, I encounter the following error:
```
SPLINE_INFO| Generating 1 splines for NONBONDED interactions
Due to 1 different atomic kinds

*******************************************************************************
* ___ *
* / \ *
* [ABORT] SPLINE_INFO| Number of points: 2346797 obtained accuracy 43.1289 *
* \___/ . MM SPLINE: no convergence on required accuracy (adjust *
* | EPS_SPLINE and rerun) *
* O/| *
* /| | *
* / \ pair_potential.F:426 *
*******************************************************************************

===== Routine Calling Stack =====

7 spline_nonbond_control
6 force_field_pack_splines
5 force_field_pack
4 force_field_control
3 fist_init
2 fist_create_force_env
1 CP2K
```
I suspect this is happening because I did not explicitly specify that each 5-atom unit is a single molecule, so CP2K is applying the Lennard-Jones potential between atoms that should be bonded. Since the atoms within each molecule are very close together, CP2K treats them as interacting via nonbonded forces, which leads to numerical instability in the Lennard-Jones calculations.

**If my assumption is correc**, How can I properly define a molecule in CP2K using the 5 carbon atoms?
How can I exclude Lennard-Jones interactions **within the same molecule**, while keeping them for interactions between different molecules?

N.B:
I understand that the carbon atoms are very close together, but this is a minimum working example meant to illustrate the issue.
I cannot change the atomic positions, as the actual problem requires keeping them fixed. Also I want to ignore the bond vibrations, rotations, torsions etc.

Any help or guidance would be greatly appreciated!

Marcella Iannuzzi

unread,

Feb 18, 2025, 9:34:48 AM2/18/25

to cp2k

Hi

If that is the problem , a straightforward way would be to pass the connectivity information through a connectivity file.

The connectivity generated by CP2K can also be tuned from

FORCE_EVAL

SUBSYS

TOPOLOGY

GENERATE

Coordinates passed in PDB format, where molecule name and molecule number are specified, might also work.

Regards

Marcella

Rajgowrav Cheenikundil

unread,

Feb 19, 2025, 12:32:22 PM2/19/25

to cp...@googlegroups.com

Hi Marcella,

Thank you for your quick response. Meanwhile I implemented your advice and replaced my "xyz" type coordinate file with a "pdb" file, with CONNEC information. But apparently CP2K is unable to read the connectivity information from there. So I defined a "PSF" file for the molecule myself and tried reading the connectivity information from there. However I'm getting the error:+
``` PSF_INFO| CMA :: Unimplemented keyword in CP2K PSF/UPSF format!```
Here is my PDB and PSF files. Could you please help me understand what is going wrong?
Is there any alternate ways to define the connectivity for artificial molecules like these? I appreciate your time and help.

1. C5_20.pdb
```
HEADER
TITLE Built with Packmol
REMARK Packmol generated pdb file
REMARK Home-Page: http://m3g.iqm.unicamp.br/packmol
REMARK
ATOM 1 C1 CC AA 1 3.763 1.559 -0.001 1.00 0.00 C
ATOM 2 C2 CC AA 1 4.432 1.969 0.619 1.00 0.00 C
ATOM 3 C3 CC AA 1 5.101 2.380 1.238 1.00 0.00 C
ATOM 4 C4 CC AA 1 5.770 2.791 1.858 1.00 0.00 C
ATOM 5 C5 CC AA 1 6.439 3.202 2.477 1.00 0.00 C
ATOM 6 C1 CC AA 2 7.024 10.035 5.233 1.00 0.00 C
ATOM 7 C2 CC AA 2 6.590 9.963 4.335 1.00 0.00 C
ATOM 8 C3 CC AA 2 6.156 9.890 3.437 1.00 0.00 C
ATOM 9 C4 CC AA 2 5.722 9.817 2.539 1.00 0.00 C
ATOM 10 C5 CC AA 2 5.288 9.744 1.641 1.00 0.00 C
ATOM 11 C1 CC AA 3 10.002 6.551 9.985 1.00 0.00 C
ATOM 12 C2 CC AA 3 9.995 5.677 9.500 1.00 0.00 C
ATOM 13 C3 CC AA 3 9.989 4.802 9.014 1.00 0.00 C
ATOM 14 C4 CC AA 3 9.983 3.928 8.529 1.00 0.00 C
ATOM 15 C5 CC AA 3 9.977 3.053 8.044 1.00 0.00 C
ATOM 16 C1 CC AA 4 3.604 2.706 9.741 1.00 0.00 C
ATOM 17 C2 CC AA 4 2.835 3.343 9.802 1.00 0.00 C
ATOM 18 C3 CC AA 4 2.067 3.980 9.863 1.00 0.00 C
ATOM 19 C4 CC AA 4 1.298 4.616 9.925 1.00 0.00 C
ATOM 20 C5 CC AA 4 0.529 5.253 9.986 1.00 0.00 C
CONECT 1 2
CONECT 2 1 3
CONECT 3 2 4
CONECT 4 3 5
CONECT 5 4
CONECT 6 7
CONECT 7 6 8
CONECT 8 7 9
CONECT 9 8 10
CONECT 10 9
CONECT 11 12
CONECT 12 11 13
CONECT 13 12 14
CONECT 14 13 15
CONECT 15 14
CONECT 16 17
CONECT 17 16 18
CONECT 18 17 19
CONECT 19 18 20
CONECT 20 19
END
```
2. C5_molecule.psf
```
PSF CMAP

1 !NTITLE

5 !NATOM
1 MOL 1 CC C1 C 0.000000 12.0110 0
2 MOL 1 CC C2 C 0.000000 12.0110 0
3 MOL 1 CC C3 C 0.000000 12.0110 0
4 MOL 1 CC C4 C 0.000000 12.0110 0
5 MOL 1 CC C5 C 0.000000 12.0110 0

4 !NBOND: bonds
1 2 2 3 3 4 4 5

3 !NTHETA: angles
1 2 3 2 3 4 3 4 5

2 !NPHI: dihedrals
1 2 3 4 2 3 4 5

0 !NIMPHI: impropers

0 !NCRTERM: cross-terms
```
3. Input file

&BOND
ATOMS C C
K 900.0
KIND HARMONIC
R0 1.0
&END BOND

&BEND
ATOMS C C C
KIND HARMONIC
K [rad^-2kcalmol] 1100.0
THETA0 [deg] 180
&END BEND

&CHARGE
ATOM C
CHARGE 0
&END CHARGE

&NONBONDED

&LENNARD-JONES
atoms C C

EPSILON 0.08
SIGMA 3.1507
RCUT 11.4

RMIN 1.0
&END LENNARD-JONES

&END NONBONDED

&END FORCEFIELD

&POISSON
PERIODIC NONE
&EWALD
EWALD_TYPE NONE
&END EWALD
&END POISSON

&END MM

&SUBSYS

&CELL
PERIODIC XYZ

ABC 30 30 30
&END CELL

&TOPOLOGY
COORD_FILE_FORMAT pdb
COORD_FILE_NAME C5_20.pdb
CONN_FILE_FORMAT PSF
CONN_FILE_NAME C5_molecule.psf

--
You received this message because you are subscribed to a topic in the Google Groups "cp2k" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cp2k/V5RlsZPw1Lk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cp2k+uns...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/cp2k/c8965abb-ad8a-41c7-98b1-ca8674207b75n%40googlegroups.com.

Reply all

Reply to author

Forward