I am wanting to simulate polymer melts. I created the initial files using the chain.f tool built in Lammps, and running the simulation protocol from the paper attached.
When equilibrated my melt Radial Dist Function matches the one given in the paper, also the end-to-end distance Re^2 parameter is matching as well, but the Gyration radius doesn't match the expected value of Rg^2 = Re^2/6, but is 10 times smaller. Regardless of the chain length, does anyone have any (even intuitive suggestion) what may have gone wrong?
http://lammps.sandia.gov/threads/pdfEfEwE2ix6N.pdf
How do you calculated gyration radius? Is it any built-in procedure? Do you particles have any effective "mass"? Was it taken into account calculating R_g?
Thank you all for the replies.
When equilibrated my melt exhibits correct end-to-end distance (matching the theoretical = C * l^2 N_monomers ( I calculate this parameter using VMD, MDanalysis tools)
It also matched the RDF provided in the paper.
But the ratio of / which is supposed to be equal to 6 is equal to 0.6 ( I calculate the using Mdanalysis, Lammps and VMD)
So my Gyration radius squared is 10 times larger then the expected value.
When I transform my data file into Real units and run a normal Lammps simulation with a forcefield with angles,dihedrals and etc. I again get the ratio which is around 0.6
It looks like its a systematic error, but I double checked the units and everything seems to be right. This issue is universal for chains of different lengths.
Here is my equilibration script:
# Kremer-Grest model.
units lj
atom_style bond
special_bonds lj/coul 0 1 1
read_data $input_data_file
neighbor 0.4 bin
neigh_modify every 1 delay 1
comm_modify vel yes
bond_style fene
bond_coeff * 30.0 1.5 1.0 1.0
dump mydump all dcd 5000 temporary_traj.dcd
timestep 0.01
thermo 100
pair_style dpd 1.0 1.0 122347
pair_coeff * * 25 4.5 1.0
velocity all create 1.0 17786140
fix 1 all nve/limit 0.05
run 500
unfix 1
fix 1 all nve
run 50000
pair_coeff * * 50.0 4.5 1.0
velocity all create 1.0 15086120
run 50
pair_coeff * * 100.0 4.5 1.0
velocity all create 1.0 15786120
run 50
pair_coeff * * 150.0 4.5 1.0
velocity all create 1.0 15486120
run 50
pair_coeff * * 200.0 4.5 1.0
velocity all create 1.0 17986120
run 100
pair_coeff * * 250.0 4.5 1.0
velocity all create 1.0 15006120
run 100
pair_coeff * * 500.0 4.5 1.0
velocity all create 1.0 15087720
run 100
pair_coeff * * 1000.0 4.5 1.0
velocity all create 1.0 15086189
run 100
pair_style hybrid/overlay lj/cut 1.122462 dpd/tstat 1.0 1.0 1.122462 122347
pair_modify shift yes
pair_coeff * * lj/cut 1.0 1.0 1.122462
pair_coeff * * dpd/tstat 4.5 1.122462
velocity all create 1.0 1508612013
run 50
velocity all create 1.0 15086121
run 50
velocity all create 1.0 15086111
run 50
velocity all create 1.0 15086125
run 1000000
I also created a small molecule ( 8 atoms) and manually tested the calculation, the scripts are giving me correct values. My mass for LJ units is 1.0, and for the Real ones 14.0
Thank you Dr. Alavi,
Yeah indeed you are right, for some reason my gyration was the gyration radius of the whole system, whereas the R2 was the property of a single molecule.
Attached are snapshotst of my melt
Hi Xinmeng,
Actually I don't know where to find it, I made a mistake calling it open source, since its just free to use. What do you need it for? For me since I didn't need to modify it, i didn't need to have the source.
For programms where you can see the source I would recommend to use Lammps/tools/ : moltemplate, chain.f, polybond.
Moltemplate has a capability to produce structures of any difficulty.
Vasiliy
- Badges
- Science topic
- Similar topics
- Chemistry
- Theoretical Chemistry
- Computational Chemistry
- Molecular Dynamics