- Posts: 30
cg topology
- msaeedi
-
Topic Author
- Offline
- Junior Boarder
Less
More
10 years 6 months ago #2320
by msaeedi
cg topology was created by msaeedi
Hi,
I mapped a molecule with 3 cg beads and determined bonded parameters with compared by AA simulation. After that I wanted to measure partitioning free energy for AA and CG models. I determined Free energy in water with g_bar method .These values are -1.2 and 57.5 KJ/mol for AA and CG respectively. These values are different from each other. I changed cg bead types and even though I changed the mapping scheme but the free energy in cg model don't change very much. What is wrong in my work? What can I do to solve this?
I mapped a molecule with 3 cg beads and determined bonded parameters with compared by AA simulation. After that I wanted to measure partitioning free energy for AA and CG models. I determined Free energy in water with g_bar method .These values are -1.2 and 57.5 KJ/mol for AA and CG respectively. These values are different from each other. I changed cg bead types and even though I changed the mapping scheme but the free energy in cg model don't change very much. What is wrong in my work? What can I do to solve this?
Please Log in or Create an account to join the conversation.
- xavier
-
- Offline
- Admin
Less
More
- Posts: 416
10 years 6 months ago #2518
by xavier
Replied by xavier on topic cg topology
It is a bit difficult to help you with s little information. Many things might get wrong at may stages ... may be your CG run is wrong or the topology wrong, but may be your atomistic data is wrong ...
You need to give more detail on your system and protocol to get more feedback.
You need to give more detail on your system and protocol to get more feedback.
Please Log in or Create an account to join the conversation.
- msaeedi
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 30
10 years 6 months ago - 10 years 6 months ago #2520
by msaeedi
Replied by msaeedi on topic cg topology
Hi,
I worked with fluorouracil molecule (C4H3FN2O2) that it is a pyrimidine analog. I guessed that it had 3 beads with name P3, P3 and C5 that correspond to NHCO, NHCO and CH2CHF groups. Then, I simulated fluorouracil with 2167 spc216 water model with gromos43al force field by EM/NVT/NPT/MD process. I simulated this system for 100ns. After that, I used g_dist to determine bond distance for each 2 groups (NHCO,NHCO),(NHCO,CH2CHF). Then, I inserted this bond distance in initial cg topology. This topology is:
;;;;;; FLU
[moleculetype]
; molname nrexcl
FLU 1
[atoms]
; id type resnr residu atom cgnr charge
1 P3 1 FLU C1 1 0
2 P3 1 FLU C2 2 0
3 C5 1 FLU C3 3 0
[bonds]
; i j funct length force.c.
1 2 1 0.284 20000
1 3 1 0.248 20000
2 3 1 0.336 20000
I used molmaker.py to transfer this to gro file.
Then I simulated cg molecule with 545 cg waters. I did simulations in EM/NPT/MD process and I simulated this molecule for 100 ns. After that, I measured bond distance with g_dist and found corresponding AA bond distance. I thought that this cg topology was right.
For verifying cg topology, I wanted to measure partitioning free energy. I used g_bar method to determine free energy of salvation in water. I simulated aa system with EM(steep)/NVT(sd)/NPT(sd)/MD(sd) process for each lambda. After that, I simulated cg system with EM(steep)/NPT(md)/MD(md) process for each lambda and then I determined free energy. They were -1.2 KJ/Mol for AA and 57.7 KJ/mol for CG model. These values are very different from each other. May I worked wrong that got these values?
I worked with fluorouracil molecule (C4H3FN2O2) that it is a pyrimidine analog. I guessed that it had 3 beads with name P3, P3 and C5 that correspond to NHCO, NHCO and CH2CHF groups. Then, I simulated fluorouracil with 2167 spc216 water model with gromos43al force field by EM/NVT/NPT/MD process. I simulated this system for 100ns. After that, I used g_dist to determine bond distance for each 2 groups (NHCO,NHCO),(NHCO,CH2CHF). Then, I inserted this bond distance in initial cg topology. This topology is:
;;;;;; FLU
[moleculetype]
; molname nrexcl
FLU 1
[atoms]
; id type resnr residu atom cgnr charge
1 P3 1 FLU C1 1 0
2 P3 1 FLU C2 2 0
3 C5 1 FLU C3 3 0
[bonds]
; i j funct length force.c.
1 2 1 0.284 20000
1 3 1 0.248 20000
2 3 1 0.336 20000
I used molmaker.py to transfer this to gro file.
Then I simulated cg molecule with 545 cg waters. I did simulations in EM/NPT/MD process and I simulated this molecule for 100 ns. After that, I measured bond distance with g_dist and found corresponding AA bond distance. I thought that this cg topology was right.
For verifying cg topology, I wanted to measure partitioning free energy. I used g_bar method to determine free energy of salvation in water. I simulated aa system with EM(steep)/NVT(sd)/NPT(sd)/MD(sd) process for each lambda. After that, I simulated cg system with EM(steep)/NPT(md)/MD(md) process for each lambda and then I determined free energy. They were -1.2 KJ/Mol for AA and 57.7 KJ/mol for CG model. These values are very different from each other. May I worked wrong that got these values?
Last edit: 10 years 6 months ago by msaeedi.
Please Log in or Create an account to join the conversation.
- xavier
-
- Offline
- Admin
Less
More
- Posts: 416
10 years 6 months ago #2536
by xavier
All the procedure looks fine to me but I agree with you that the results are not very convincing!
I assume you did not consider all the EM/NVT/NPT/MD data in the g_bar calculation but only the part of MD that would show convergence. Did you check convergence at each lambda point?
It does not make much sense to have such a high value in the CG calculation. The beads you selected are mostly hydrophobic ... I would check these runs first.
Which input files did you use to generate your data? You need to use a soft core potential ... did you?
Replied by xavier on topic cg topology
msaeedi wrote: Hi,
I worked with fluorouracil molecule (C4H3FN2O2) that it is a pyrimidine analog. I guessed that it had 3 beads with name P3, P3 and C5 that correspond to NHCO, NHCO and CH2CHF groups. Then, I simulated fluorouracil with 2167 spc216 water model with gromos43al force field by EM/NVT/NPT/MD process. I simulated this system for 100ns. After that, I used g_dist to determine bond distance for each 2 groups (NHCO,NHCO),(NHCO,CH2CHF). Then, I inserted this bond distance in initial cg topology. This topology is:
;;;;;; FLU
[moleculetype]
; molname nrexcl
FLU 1
[atoms]
; id type resnr residu atom cgnr charge
1 P3 1 FLU C1 1 0
2 P3 1 FLU C2 2 0
3 C5 1 FLU C3 3 0
[bonds]
; i j funct length force.c.
1 2 1 0.284 20000
1 3 1 0.248 20000
2 3 1 0.336 20000
I used molmaker.py to transfer this to gro file.
Then I simulated cg molecule with 545 cg waters. I did simulations in EM/NPT/MD process and I simulated this molecule for 100 ns. After that, I measured bond distance with g_dist and found corresponding AA bond distance. I thought that this cg topology was right.
For verifying cg topology, I wanted to measure partitioning free energy. I used g_bar method to determine free energy of salvation in water. I simulated aa system with EM(steep)/NVT(sd)/NPT(sd)/MD(sd) process for each lambda. After that, I simulated cg system with EM(steep)/NPT(md)/MD(md) process for each lambda and then I determined free energy. They were -1.2 KJ/Mol for AA and 57.7 KJ/mol for CG model. These values are very different from each other. May I worked wrong that got these values?
All the procedure looks fine to me but I agree with you that the results are not very convincing!
I assume you did not consider all the EM/NVT/NPT/MD data in the g_bar calculation but only the part of MD that would show convergence. Did you check convergence at each lambda point?
It does not make much sense to have such a high value in the CG calculation. The beads you selected are mostly hydrophobic ... I would check these runs first.
Which input files did you use to generate your data? You need to use a soft core potential ... did you?
Please Log in or Create an account to join the conversation.
- msaeedi
-
Topic Author
- Offline
- Junior Boarder
Less
More
- Posts: 30
10 years 5 months ago #2612
by msaeedi
Replied by msaeedi on topic cg topology
Hi. Thanks very much for your response and concentration.
When I wanted to calculate free energy calculation with g_bar ,I used em.gro for NVT and nvt.gro for NPT and npt.gro for MD simulation of each lambda point. I checked convergences of energy for each lambda.
For example I used this md.mdp file for aa simulation. I copied this file in Gromacs tutorial of " calculation of free energy for methane in water".
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 1000000 ; 2 ns
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxtcout = 0
xtc-precision = 1000
; Neighborsearching and short-range nonbonded interactions
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.0
; Electrostatics
coulombtype = PME
rcoulomb = 1.0
; van der Waals
vdw-type = switch
rvdw-switch = 0.8
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
optimize_fft = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 300
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 0.5
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda = 0.0
delta_lambda = 0
foreign_lambda = 0.05
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
couple-moltype = FLU ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; Do not generate velocities
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
for CG simulation, I copied free energy section in mdp file. This md.mdp file in cg simulation for lamnda=0 is:
; VARIOUS PREPROCESSING OPTIONS =
title = Martini
cpp = /usr/bin/cpp
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0.0
dt = 0.020
nsteps = 500000
; number of steps for center of mass motion removal =
nstcomm = 1
comm-grps =
; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 5000
nstvout = 5000
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 1000
nstenergy = 100
; Output frequency and precision for xtc file =
nstxtcout = 1000
xtc_precision = 100
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps =
; Selection of energy groups =
energygrps = FLU W
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 10
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or none =
pbc = xyz
; nblist cut-off =
rlist = 1.2
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Shift
rcoulomb_switch = 0.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
; Method for doing Van der Waals =
vdw_type = Shift
; cut-off lengths =
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = Berendsen
; Groups to couple separately =
tc-grps = FLU W
; Time constant (ps) and reference temperature (K) =
tau_t = 1.0 1.0
ref_t = 310 310
; Pressure coupling =
Pcoupl = berendsen
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p = 0.2
compressibility = 3e-5
ref_p = 1.0
; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel = no
gen_temp = 105
gen_seed = 473529
; Free energy control stuff
free_energy = yes
init_lambda = 0.0
delta_lambda = 0
foreign_lambda = 0.05
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
couple-moltype = FLU ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; OPTIONS FOR BONDS =
constraints = none
; Type of constraint algorithm =
constraint_algorithm = Lincs
; Do not constrain the start configuration =
unconstrained_start = no
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle = 30
also, I used initial gro file "FLU in water" for simulation of each lambda point. I think that I don’t use right coulomb and vdwd type of potential in cg simulation for free energy calculation. If it is possible for you, please help me.
When I wanted to calculate free energy calculation with g_bar ,I used em.gro for NVT and nvt.gro for NPT and npt.gro for MD simulation of each lambda point. I checked convergences of energy for each lambda.
For example I used this md.mdp file for aa simulation. I copied this file in Gromacs tutorial of " calculation of free energy for methane in water".
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 1000000 ; 2 ns
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxtcout = 0
xtc-precision = 1000
; Neighborsearching and short-range nonbonded interactions
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.0
; Electrostatics
coulombtype = PME
rcoulomb = 1.0
; van der Waals
vdw-type = switch
rvdw-switch = 0.8
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
optimize_fft = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 300
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 0.5
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda = 0.0
delta_lambda = 0
foreign_lambda = 0.05
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
couple-moltype = FLU ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; Do not generate velocities
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
for CG simulation, I copied free energy section in mdp file. This md.mdp file in cg simulation for lamnda=0 is:
; VARIOUS PREPROCESSING OPTIONS =
title = Martini
cpp = /usr/bin/cpp
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0.0
dt = 0.020
nsteps = 500000
; number of steps for center of mass motion removal =
nstcomm = 1
comm-grps =
; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 5000
nstvout = 5000
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 1000
nstenergy = 100
; Output frequency and precision for xtc file =
nstxtcout = 1000
xtc_precision = 100
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps =
; Selection of energy groups =
energygrps = FLU W
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 10
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or none =
pbc = xyz
; nblist cut-off =
rlist = 1.2
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Shift
rcoulomb_switch = 0.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
; Method for doing Van der Waals =
vdw_type = Shift
; cut-off lengths =
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = Berendsen
; Groups to couple separately =
tc-grps = FLU W
; Time constant (ps) and reference temperature (K) =
tau_t = 1.0 1.0
ref_t = 310 310
; Pressure coupling =
Pcoupl = berendsen
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p = 0.2
compressibility = 3e-5
ref_p = 1.0
; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel = no
gen_temp = 105
gen_seed = 473529
; Free energy control stuff
free_energy = yes
init_lambda = 0.0
delta_lambda = 0
foreign_lambda = 0.05
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
couple-moltype = FLU ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
nstdhdl = 10
; OPTIONS FOR BONDS =
constraints = none
; Type of constraint algorithm =
constraint_algorithm = Lincs
; Do not constrain the start configuration =
unconstrained_start = no
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle = 30
also, I used initial gro file "FLU in water" for simulation of each lambda point. I think that I don’t use right coulomb and vdwd type of potential in cg simulation for free energy calculation. If it is possible for you, please help me.
Please Log in or Create an account to join the conversation.
Time to create page: 0.163 seconds