Colvars / colvars

Collective variables library for molecular simulation and analysis programs
http://colvars.github.io/
GNU Lesser General Public License v3.0
201 stars 59 forks source link

Why is my ABF simulation sampling from a small set of values? #592

Closed wisecashew closed 10 months ago

wisecashew commented 11 months ago

I am running a simulation of a polymer in water. I am trying to calculate the free energy of the radius of gyration using ABF. I am using LAMMPS. I would like an F(Rg) curve as 5<=Rg<=30, but the ends can be flexible. This is my input file for the simulation:

Colvarstrajfrequency    1000
Colvarsrestartfrequency 20000

colvar {
    name my_rg
    gyration {
        atoms { atomNumbersRange 1-572 }
    }
    lowerBoundary 5.0 
    upperBoundary 30.0
    # width=0.5
}

harmonicWalls {
    name              wall
    colvars           my_rg
    lowerWalls        5.0 
    upperWalls        30.0
    lowerWallConstant 25.0
    upperWallConstant 25.0
}

abf {
    name abf_r
    colvars     my_rg
    outputFreq  5000
    fullSamples 1000
    historyFreq 50000
}

On running the simulation, I looked at out.pmf file and I see:

# 1
#  4.50000000000000e+00  1.00000000000000e+00         26  0

  5.00000000000000e+00   1.20005097266126e+02
  6.00000000000000e+00   1.20005097266126e+02
  7.00000000000000e+00   1.20005097266126e+02
  8.00000000000000e+00   1.20005097266126e+02
  9.00000000000000e+00   2.81986236616591e+01
  1.00000000000000e+01   0.00000000000000e+00
  1.10000000000000e+01   0.00000000000000e+00
  1.20000000000000e+01   0.00000000000000e+00
  1.30000000000000e+01   0.00000000000000e+00
  1.40000000000000e+01   0.00000000000000e+00
  1.50000000000000e+01   0.00000000000000e+00
  1.60000000000000e+01   0.00000000000000e+00
  1.70000000000000e+01   0.00000000000000e+00
  1.80000000000000e+01   0.00000000000000e+00
  1.90000000000000e+01   0.00000000000000e+00
  2.00000000000000e+01   0.00000000000000e+00
  2.10000000000000e+01   0.00000000000000e+00
  2.20000000000000e+01   0.00000000000000e+00
  2.30000000000000e+01   0.00000000000000e+00
  2.40000000000000e+01   0.00000000000000e+00
  2.50000000000000e+01   0.00000000000000e+00
  2.60000000000000e+01   0.00000000000000e+00
  2.70000000000000e+01   0.00000000000000e+00
  2.80000000000000e+01   0.00000000000000e+00
  2.90000000000000e+01   0.00000000000000e+00
  3.00000000000000e+01   0.00000000000000e+00

It seems like only integer values less than 9 have been sampled. Furthermore, if I look at out.count:

# 1
#  5.00000000000000e+00  1.00000000000000e+00         25  0

  5.50000000000000e+00                      0   
  6.50000000000000e+00                      0   
  7.50000000000000e+00                      0   
  8.50000000000000e+00               12040750
  9.50000000000000e+00                7959250
  1.05000000000000e+01                      0   
  1.15000000000000e+01                      0   
  1.25000000000000e+01                      0   
  1.35000000000000e+01                      0   
  1.45000000000000e+01                      0   
  1.55000000000000e+01                      0   
  1.65000000000000e+01                      0   
  1.75000000000000e+01                      0   
  1.85000000000000e+01                      0   
  1.95000000000000e+01                      0   
  2.05000000000000e+01                      0   
  2.15000000000000e+01                      0   
  2.25000000000000e+01                      0   
  2.35000000000000e+01                      0   
  2.45000000000000e+01                      0   
  2.55000000000000e+01                      0   
  2.65000000000000e+01                      0   
  2.75000000000000e+01                      0   
  2.85000000000000e+01                      0   
  2.95000000000000e+01                      0   

It seems to confirm that only about Rg from 8-9 are being sampled. I decided to restart the run using the previous state file and this input file:

Colvarstrajfrequency    1000
Colvarsrestartfrequency 20000

colvar {
    name my_rg
    gyration {
        atoms { atomNumbersRange 1-572 }
    }
    lowerBoundary 5.0 
    upperBoundary 30.0
    # width=0.5
}

harmonicWalls {
    name              wall
    colvars           my_rg
    lowerWalls        5.0 
    upperWalls        30.0
    lowerWallConstant 25.0
    upperWallConstant 25.0
}

abf {
    name abf_r
    # inputPrefix out 
    colvars     my_rg
    outputFreq  5000
    fullSamples 1000
    historyFreq 50000
}

But I see no improvements in the values that are being sampled. How do I improve the sampling from the simulation? Are there any ways to tweak the force parameters for better sampling? I am running metadynamics simulations using PLUMED to check that both Free Energy Surfaces match, and as is it stands, it seems like the ABF is simply not sampling well enough.

I would greatly appreciate any advice you have for me!

jhenin commented 11 months ago

I think the grid on which you discretize the free energy is too coarse, so the ABF bias cannot do its job. You need to make the "width" parameter of your collective variable smaller.

giacomofiorin commented 10 months ago

@wisecashew I hope that changing the width parameter was useful. Although it's unlikely that ABF and metadynamics would yield significantly different results, do keep in mind that there is a metadynamics implementation for Colvars as well: https://colvars.github.io/colvars-refman-lammps/colvars-refman-lammps.html#sec:colvarbias_meta