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

[Planck]

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 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/29a20e23-538c-4d81-95bc-cb6020882f5an%40googlegroups.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20250218/568f1cd5/attachment.htm>