[CP2K-user] [CP2K:21174] Re: CP2K: How to Define a Molecule and Restrict Lennard-Jones Interactions Within It?

Rajgowrav Cheenikundil rajg.chnk at gmail.com
Wed Feb 19 17:31:24 UTC 2025


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

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


      &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


    &END TOPOLOGY

  &END SUBSYS

  STRESS_TENSOR NUMERICAL
&END FORCE_EVAL
```






On Tue, Feb 18, 2025 at 4:34 PM Marcella Iannuzzi <marci.akira at gmail.com>
wrote:

> 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 <https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL.html>
>
>
>
>    - SUBSYS
>    <https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS.html>
>
>
>
>    - TOPOLOGY
>    <https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/TOPOLOGY.html>
>
>
>
>    - GENERATE
>
>
> Coordinates passed in PDB format, where molecule name and molecule number
> are specified, might also work.
>
> Regards
> Marcella
>
> On Tuesday, February 18, 2025 at 3:06:02 PM UTC+1 Planck wrote:
>
>> 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!
>>
>> --
> 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+unsubscribe at googlegroups.com.
> To view this discussion visit
> https://groups.google.com/d/msgid/cp2k/c8965abb-ad8a-41c7-98b1-ca8674207b75n%40googlegroups.com
> <https://groups.google.com/d/msgid/cp2k/c8965abb-ad8a-41c7-98b1-ca8674207b75n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>

-- 
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+unsubscribe at googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/cp2k/CAFgiMy-_FJd1Yuhg-LF%3DruX%2BRfQ-4zXpvQdZna0-yNYxpF-qsQ%40mail.gmail.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20250219/e5d8c378/attachment-0001.htm>


More information about the CP2K-user mailing list