Montecarlo

Matt McGrath obfis... at gmail.com
Wed Nov 10 07:14:45 UTC 2010


Hello Francesco.  Let me try to understand exactly what you're looking
for, and let's see if CP2K can do it.

First, have you looked at the MC section in input_cp2k_motion.F?  If
you search for the MC keywords in that file, you can find a short
description of them.  For example, PMVOLUME gives "The probability of
attempting a volume move."  In your original input file, this is equal
to 1.0.  Move probabilities are calculated somewhat obtusely, due to
legacy reasons with another code (you can look in mc_ensembles.F for
the details...it's not difficult, just odd), but because this is equal
to 1.0 and PMSWAP, PMTRAION, and PMTRANS are all equal to zero , you
are doing 100% volume displacements.  All PM variables are some type
of probability.  If you want to do just conformational changes,
PMTRAION should be equal to 1.0, and the other three should be equal
to 0.0.

Now, all this does is perform a Monte Carlo simulation using only
conformation changes (canonical ensemble).  That is, it selects a
degree of freedom at random, attempts to move it by a random amount
(between 0 and RMBOND/RMBEND/RMDIHEDRAL, depending on which degree of
freedom selected), calculates the acceptance based on the standard
Metropolis rule, and either accepts or rejects.  That's it.  If that's
exactly what you're looking for, then great.  If you're looking for
something else, CP2K-MC probably can't do it right now.  What,
exactly, do you mean by "conformational analysis"?

No, AVBMC won't help you for a single molecule.  AVBMC is most often
used in vapor phases for associating fluids which tend to cluster
(water, alcohols, HF, acetic acid...).  It's a method to cross the
energy barrier in creating and destroying a cluster.

                           Cheers, Matt


On Nov 9, 5:09 pm, Francesco <francesc... at gmail.com> wrote:
> HiMatt,
>
> Thanks for your psf file, but I cannot obtain the conformational
> analysis for a single ethane molecule, with any attempts I have done.
>
> After 1000 steps I obtain just small volume changes, that of course,
> for a single molecule is not useful.
> If I block the volume variations, in the output file ethane.inp.moves
> I obtain just empty_avbmc movements, i.e. the program does nothing.
>
> I started to look in the literature for some information about AVBMC
> method: is it suitable for conformational analysis also for single
> molecules? Every application I found is about condensed phase.
>
> Could you paste the whole input file? I guess I cannot understand how
> to use some keyword.
>
> Thank you in advance,
> Francesco
>
> On Sep 15, 9:53 am,Matt<obfis... at gmail.com> wrote:
>
>
>
> > Hello Francesco.
>
> > One of the quirks with MC in CP2K is that it needs topology files
> > defined for every molecule type.  The reason is because we noticed
> > that the code would crash sometimes during a regular NVT run when CP2K
> > found that the number of molecules had changed (the automatic
> > molecular determination).
>
> > I have everything working with the following .psf file and topology
> > section:
>
> >     &TOPOLOGY
> >       CONNECTIVITY MOL_SET
> >       &MOL_SET
> >         &MOLECULE
> >           NMOL 1
> >           CONN_FILE_NAME topology_atoms_ETH.psf
> >         &END
> >       &END
> >     &END
>
> > Make sure that the coordinates in the input file match the order in
> > the topology file (i.e. hydrogens 3,4,5 are bonded to carbon 1).  This
> > is a rather clumsy way to do it, but it's the best I could come up
> > with when designing the routines.  The most difficult part of this is
> > creating the psf file: the formatting is VERY strict (number of
> > spaces), and a slight mistake will cause a seg fault (as I discovered
> > when creating this file).  I took the format from cp2k/tests/Fist/
> > sample_psf/solute1.psf.  Good luck!
>
> >                                   Cheers,Matt
>
> >  PSF EXT
>
> >          1 !NTITLE
> >      Topology file for ethane
>
> >          8 !NATOM
> >          1 ETH             1 ETH      C       C       0.000000
> > 12.000000        0
> >          2 ETH             1 ETH      C       C       0.000000
> > 12.000000        0
> >          3 ETH             1 ETH      H       H       0.000000
> > 1.008000        0
> >          4 ETH             1 ETH      H       H       0.000000
> > 1.008000        0
> >          5 ETH             1 ETH      H       H       0.000000
> > 1.008000        0
> >          6 ETH             1 ETH      H       H       0.000000
> > 1.008000        0
> >          7 ETH             1 ETH      H       H       0.000000
> > 1.008000        0
> >          8 ETH             1 ETH      H       H       0.000000
> > 1.008000        0
>
> >          7 !NBOND
> >          1         2         1         3         1         4
> > 1         5
> >          2         6         2         7         2         8
>
> >          6 !NTHETA
> >          3         1         2         4         1         2
> > 5         1         2
> >          6         2         1         7         2         1
> > 8         2         1
>
> >          9 !NPHI
> >          3         1         2         6         4         1
> > 2         6
> >          5         1         2         6         3         1
> > 2         7
> >          4         1         2         7         5         1
> > 2         7
> >          3         1         2         8         4         1
> > 2         8
> >          5         1         2         8
>
> >          0 !NIMPHI
>
> >          0 !NDON
>
> >          0 !NACC
>
> >          0 !NNB
>
> >          0 !NGRP
>
> > On Sep 14, 7:38 pm, Francesco <francesc... at gmail.com> wrote:
>
> > > Hi to all,
>
> > > I would like to use the CP2K-Montecarlomodule in order to do a
> > > conformational analysis test, with DFT.
>
> > > In order to understand what to do, I have started with ethane (and
> > > other small molecules), but any attempt failed till now. I've used the
> > > file in the regtest as a template, but in the best case I obtained the
> > > error attached below.
>
> > > After some tests, it seems to me that the Number of molecule types
> > > vary as the number of atoms: I tried also to define a psf file, but
> > > with other errors. By the way with "CONN_FILE_FORMAT OFF" should
> > > overcome this request. Or am I wrong?
>
> > > Has anybody of you done something similar?
> > > Any help will be appreciated.
>
> > > Thanks in advance,
> > > Francesco
>
> > > --------------------------------------------------------------------------- -----
> > >  Number of molecule types            8
> > >  &MOLECULE sections found            1
>
> > >  *
> > >  *** ERROR in mc_types.F/mc_input_file_create  ***
> > >  *
>
> > >  *** Did not find MOLECULE sections for every molecule in the
> > > simulation  ***
> > >  *** (make sure both input files have all
> > > types)                          ***
>
> > > --------------------------------------------------------------------------- -------
> > > &FORCE_EVAL
> > >   METHOD Quickstep
> > >   &DFT
> > >     &MGRID
> > >       CUTOFF 100
> > >     &END MGRID
> > >     &QS
> > >       EXTRAPOLATION USE_PREV_WF
> > >     &END QS
> > >     &SCF
> > >       SCF_GUESS ATOMIC
> > >     &END SCF
> > >     &XC
> > >       &XC_FUNCTIONAL Pade
> > >       &END XC_FUNCTIONAL
> > >       &XC_GRID
> > >         XC_DERIV SPLINE2
> > >         XC_SMOOTH_RHO NONE
> > >       &END XC_GRID
> > >     &END XC
> > >   &END DFT
> > >   &SUBSYS
> > >     &CELL
> > >       ABC 5.04977 5.04977 5.04977
> > >     &END CELL
> > >     &COORD
> > >  C     0.000000     0.000000     0.000000
> > >  C     0.000000     0.000000     1.450000
> > >  H     1.026720     0.000000     1.812996
> > >  H    -0.513360     0.889165     1.813000
> > >  H    -0.513360    -0.889165     1.813000
> > >  H    -1.026719     0.000000    -0.363000
> > >  H     0.513360     0.889165    -0.363000
> > >  H     0.513360    -0.889165    -0.363000
> > >     &END COORD
> > >     &KIND H
> > >       BASIS_SET DZVP-GTH-Pade
> > >       POTENTIAL GTH-BLYP-q1
> > >     &END KIND
> > >     &KIND C
> > >       BASIS_SET DZVP-GTH-Pade
> > >       POTENTIAL GTH-BLYP-q4
> > >     &END KIND
> > >     &TOPOLOGY
> > >       &MOL_SET
> > >         &MOLECULE
> > >           NMOL 1
> > >           CONN_FILE_FORMAT OFF
> > >         &END MOLECULE
> > >       &END MOL_SET
> > >       CONNECTIVITY MOL_SET
> > >     &END TOPOLOGY
> > >   &END SUBSYS
> > > &END FORCE_EVAL
> > > &GLOBAL
> > >   PROJECT ethane_MC
> > >   RUN_TYPE MC
> > >   PRINT_LEVEL LOW
> > > &END GLOBAL
> > > &MOTION
> > >   &MC
> > >     IUPTRANS 6400000
> > >     IUPVOLUME 3200000
> > >     LBIAS no
> > >     LSTOP yes
> > >     NMOVES 8
> > >     NSTEP 2
> > >     PMSWAP 0.0
> > >     PMTRAION 0.00
> > >     PMTRANS 0.00
> > >     PMVOLUME 1.00
> > >     PMVOL_BOX 1.00
> > >     PRESSURE 1.013
> > >     ENSEMBLE traditional
> > >     RESTART no
> > >     RESTART_FILE_NAME mc_restart_1
> > >     RMDIHEDRAL 3.0
> > >     RMANGLE 3.0
> > >     RMBOND 0.074
> > >     RMROT 26.0
> > >     RMTRANS 0.38
> > >     RMVOLUME 0.5
> > >     TEMPERATURE 398.0
> > >     AVBMC_RMIN 1.0
> > >     AVBMC_RMAX 5.0
> > >     AVBMC_ATOM 1
> > >     PBIAS 0.5
> > >     ETA 0.0
> > >     PMAVBMC_MOL 1.0
> > >     PMSWAP_MOL 1.0
> > >     PMROT_MOL 1.0
> > >     PMTRANS_MOL 1.0
> > >     PMTRAION_MOL 1.0
> > >     VIRIAL_TEMPS 300.0
> > >   &END MC
> > > &END MOTION


More information about the CP2K-user mailing list