AIMD under PBC

[Abedi, Mostafa]

Hi Everyone,

I am trying to perform AIMD for a periodic system containing 8 CH4 molecules. Everything goes well, but after a few hundred steps, methane molecules start to leave the box. This is really strange behavior because the simulation is running under PBC and if one molecule leaves the box, one of its images must enter through the opposite face with exactly the same way and direction. But this is not happening in my simulations. So, it seems PBC is not working. I believe something is wrong in my input file, but I can not figure out what it is. I greatly appreciate it if someone could help me to fix this issue. Many thanks.     

&GLOBAL
  PROJECT methane
  RUN_TYPE MD  
  WALLTIME 1000000
  PRINT_LEVEL  MEDIUM
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS            
    CHARGE 0
    MULTIPLICITY 1
    &MGRID
       CUTOFF [Ry] 400
    &END
    &QS
       METHOD GPW
       EPS_DEFAULT 1.0E-10
       EXTRAPOLATION ASPC
    &END
    &POISSON
       PERIODIC XYZ
       POISSON_SOLVER PERIODIC
    &END
    &SCF
      &PRINT
        &RESTART OFF
        &END
      &END                
      SCF_GUESS ATOMIC
      MAX_SCF 30
      EPS_SCF 1.0E-6
      &OT
        PRECONDITIONER FULL_SINGLE_INVERSE
        MINIMIZER CG
      &END OT
      &OUTER_SCF
        MAX_SCF 10
        EPS_SCF 1.0E-6
      &END
    &END SCF
    &XC
      &XC_FUNCTIONAL
       &PBE
       &END
      &END XC_FUNCTIONAL
      &VDW_POTENTIAL
         POTENTIAL_TYPE PAIR_POTENTIAL
         &PAIR_POTENTIAL
            PARAMETER_FILE_NAME dftd3.dat
            TYPE DFTD3
            REFERENCE_FUNCTIONAL PBE
            R_CUTOFF [angstrom] 16
         &END
      &END VDW_POTENTIAL
    &END XC
  &END DFT
   &SUBSYS
    &CELL
      ABC [angstrom] 9.2 9.2 9.2
      PERIODIC XYZ
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME methane.xyz
      COORD_FILE_FORMAT XYZ
    &END
    &KIND H
      ELEMENT H
      BASIS_SET SZV-MOLOPT-GTH
      POTENTIAL GTH-PBE-q1
    &END KIND
    &KIND C
      ELEMENT C
      BASIS_SET SZV-MOLOPT-GTH
      POTENTIAL GTH-PBE-q4
    &END KIND
  &END SUBSYS
  &PRINT
    &FORCES ON
     FILENAME forces
    &END
  &END
&END FORCE_EVAL

&MOTION
 &GEO_OPT
   OPTIMIZER BFGS
   MAX_ITER  1000
   MAX_DR    [bohr] 0.003
   &BFGS
   &END
 &END
  &MD
   ENSEMBLE NVE
   STEPS 10000
   TIMESTEP 0.5
   TEMPERATURE 500
 &END MD
 &PRINT
   &CELL
    FILENAME cell
   &END
 &END
&END MOTION

Best,

Mostafa

[Marcella Iannuzzi]

Dear Mostafa

CP2K does not wrap the coordinates into the primary box, but works with the absolute coordinates.

The PBC are anyway internally applied by calculating energy and forces, hence the interaction field in which your molecules move is correct. 

To verify the actual situation within the simulation box just apply PBC when using the visualisation tool of your choice.

Regards

Marcella

[Abedi, Mostafa]

Dear Marcella,

Many thanks for your help. Now, I understand it. I've found a way to wrap the coordinate into the box using VMD:

pbc set {xxx xxx xxx} -all  # the "xxx" is the length of the box

pbc box

pbc wrap -all

This works well for the NVE and NVT simulations that the length of the box is fixed, but it is not a case for the NpT simulation box, where the volume changes. Do you have any idea how that would be for the NpT one? I greatly appreciate your help.

Thanks,

Mostafa

--
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/717d7834-fa56-4517-933c-7447156d3100n%40googlegroups.com.

[Marcella Iannuzzi]

Dear Mostafa, 

The pbc tool of VMD allow to provide the information about the box per each frame.

As alternative, one can write a script to post process the cp2k output and apply the pbc as desired.

Regards

Marcella

[Krack Matthias (PSI)]

Hi Mostafa

 

You can also try the (stand-alone) program xyz2dcd (https://github.com/cp2k/cp2k/blob/master/src/motion/xyz2dcd.F) from the current cp2k trunk which dumps the frames from a CP2K xyz trajectory file in dcd format optionally with wrapped coordinates (flags –pbc and –pbc0) using either fixed cell parameters (flags –abc and –cell) or the actual cell information for each frame from the corresponding CP2K cell file of the same NpT run, e.g. by running

 

xyz2dcd –pbc –cell_file project-1.cell project-pos-1.xyz >project-pos-1_wrapped.dcd

 

The generated binary file project-pos-1.dcd in DCD format with wrapped coordinates and the corresponding cell information for each frame can be loaded into VMD or converted again back to xyz (XMOL) format using

 

dumpdcd -xyz project-pos-1.xyz –of XMOL –o project-pos-1_wrapped.xyz project-pos-1_wrapped.dcd

 

xyz2dcd and dumpdcd are compiled automatically with the current CP2K trunk version and can be found then in folder cp2k/exe/<arch>/, but they can also be compiled independently, e.g. using

 

gfortran –ffree-form –o xyz2dcd xyz2dcd.F

 

The flag –h provides further usage info.

 

HTH

 

Matthias

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/CAG%3DTvs0s4GARLBc-nf5zAq8xGsvxm3Hac%2B8ZiiGdBjwJKqJfsA%40mail.gmail.com.

[Abedi, Mostafa]

Dear Matthias,

Many thanks for your great suggestion. That is exactly what I need.

Best,

Mostafa

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/217555aed3914192a79847e146282573%40psi.ch.