[gmx-users] PMF in vacuum

[chris.neale at utoronto.ca]

Dear Eudes:
To answer your pbc vs no-pbc question, I suggest that you use pbc=no  
and set nstlist=0 rlist=0 rvdw=0 rcoulomb=0 so that you calculate all  
interactions in direct space with no cutoffs.

## Major comments that you should still investigate

1. There is no need to use a virtual atom, the pull code will use the  
center of mass. I strongly suggest that you stop using a virtual atom  
and just use the entire benzene as an argument to the pull code group.  
I have had difficulties with slightly more complicated setups of this  
type.

2. In test2.jpg, the system without pbc shows a flat PMF after the  
cut-off -- exactly what one would expect. The pbc system shows  
continued interaction -- again what I would expect. So there is  
nothing actually all that strange here. One would not expect to see  
such a drastic difference in a high dielectric such as water, but in  
vacuum I suspect that this is expected.

3. Please clarify what your cutoff was. I don't see a cut-off listed  
in your .mdp options and leaving this to the default of 1.0 nm is a  
bad idea because it can lead to confusion a times like this. I might  
assume that it was 0.7 nm based on test2.jpg, but then see the point  
#3 below.

4a. I have no idea what -DPOSRES is actually doing for you since I  
can't see your topology.

4b. Are you sure that "pull_dim = N N Y" is really what you want?  
Sometimes one wants to average over X and Y, but I am not sure that  
you do in this case.

4c. What exactly do you believe pull_r0 and pull_r1 are doing for you?

### More minor notes:

5. regarding test1.jpg: a PMF is correct to an additive constant,  
meaning that you can shift two PMFs relative to one another. These 2  
PMFs are therefore less different than they appear in your compaison  
plot, but they do differ in the slope between 1.0 - 2.0 nm. This is  
probably just a convergence issue and you will always need to do tests  
like this.

6. regarding histo.png: can you confirm that the few very short  
gaussians are due to less sampling in a few windows? In any event, the  
overlap looks good.

Chris

-- original message --

Hello Chris, thanks for your attention.
I'm sending you some links to some tests
I performed. As I said you will notice that
depending on the parameter used my simulation
shows PMF profiles quite different. Especially what
concerns to the difference between the use or not of the PBC.

https://sites.google.com/site/fileti/files/test1.jpg
https://sites.google.com/site/fileti/files/test2.jpg
https://sites.google.com/site/fileti/files/histo.png

I have constructed two very similar topologies (ben-a.itp and ben-b.itp)
where I put a virtual site in the center of benzene.
This sites were restrained to keep my molecules
fixed distance desired.

The basic details of the simulations are given below:1000

define                   = -DPOSRES
integrator               = sd
tinit                    = 0
dt                       = 0.002
nsteps                   = 5000000 or 500000
comm-mode                = angular
nstcomm                  = 1
comm-grps                = System
bd-fric                  = 1
ld-seed                  = 1993
nstlist                  = 5
ns_type                  = simple
pbc                      = no or xyz
periodic_molecules       = no
rlist                    = 0
coulombtype              = Cut-off
vdw-type                 = Cut-off
rvdw                     = 0
DispCorr                 = no
Tcoupl                   = Nose-Hoover
tc_grps                  = system
tau_t                    = 0.1
ref_t                    = 300
Pcoupl                   = no
pcoupltype               = isotropic
tau_p                    = 1.0
compressibility          = 4.5e-5
ref_p                    = 1.0
constraints              = all-bonds
constraint_algorithm     = lincs

; COM PULLING
pull                     = umbrella
pull_geometry            = distance
pull_dim                 = N N Y
pull_r1                  = 1
pull_r0                  = 1.5
pull_constr_tol          = 1e-06
pull_start               = yes
pull_nstxout             = 10
pull_nstfout             = 10
pull_ngroups             = 1
pull_group0              = BENa
pull_weights0            =
pull_pbcatom0            = 0
pull_group1              = BENb
pull_weights1            =
pull_pbcatom1            = 0
pull_vec1                = 0 0 1
pull_init1               = 0
pull_rate1               = 0.0
pull_k1                  = 1700 or