Greetings,
I am trying to simulate protein in water to determine mechanical properties. This is something new for me coming, working mostly on Crystal material. From what read, Reaxff is a good choice since it can simulate bonds dissociation. I am using lammps since it already includes the reaxff potential. The problem that usually my system consists of huge number of atoms (including water molecules). When I ran the simulation, I got an error related to the memory allocation and tge numbers of atoms in my system ( that is what I understood after searching the error online). So now I am thinking of using hybrid style although I am not sure if it will work or not but my concern is it reasonable to use reaxff forcefield for protein and another forcefield for water. I appreciate any feedback on this topic.
Regards
Hi
I have two points that I think will help you.
1. The hybrid pair style in lammps does not accept the use of 2 reaxff potentials.
2. Did you try the safezone and minicap keywords in your pair_style command ? Using these usually fixes the memory allocation problems. I had such problem even with less than 1000 atoms. Reaxff memory allocation scheme is not the same as other potentials in lammps
Hi Khaled Badawy ,
Thanks for your reply. I tried your suggestion but for some reason lammps does not recognize the 2 keywords. below i am attaching my code and the error message i received:
Input file:
units real
atom_style full
read_data data.box
pair_style reaxff lmp_control checkqeq yes safezone 1.6 minicap 100
pair_coeff * * ffield.reax C H H N O O
fix 1 all qeq/reaxff 1 0.0 10.0 1e-6 reaxff
timestep 0.0001
neighbor 1 bin
neigh_modify every 1 delay 0 check yes
thermo_style custom step temp press etotal pe vol
dump traj all xtc 50 file.xtc
fix 6 all nve
min_style cg
thermo 1
minimize 1.0e-4 1.0e-10 10000 10000
unfix 6
velocity all create 300 5812775
fix 2 all nvt temp 300.0 300.0 $(100.0*dt)
thermo 10
thermo_modify flush yes
dump traj1 all xtc 50 file1.xtc
run 30000
unfix 2
Error:
ERROR: Illegal pair_style reaxff command (src/REAXFF/pair_reaxff.cpp:265)
Last command: pair_style reaxff lmp_control checkqeq yes safezone 1.6 minicap 100
This is the error i get with the default memory allocation:
ERROR on proc 4: step 505: bondchk failed: i=509 end(i)=12751 str(i+1)=12750
Then i tried something different, when i set checkqeq no, i dont get this error anymore, however i am not sure if simulating the system without the charge update would be reasonable. Also when i include the charge equilibration, during the minimization i get too many warning about qeq not converging althtough i already prepared the structure using gromacs, energy minimization, nvt and npt so the structure should be good:
WARNING: Fix qeq/reaxff CG convergence failed after 200 iterations at step 0 (src/REAXFF/fix_qeq_reaxff.cpp:746)
So any feedback or idea will be much appreciated
By the way, there are two atom_style in lammps for reaxff, (reaxff and reax/c), I read in the lammps documentation but i couldnt understand the difference between these two.
Hi
1. reaxff and reax/c differ in the language of the implementation of the code. However, both do the same thing.
2. The qeq minimization is very tricky. I believe that using the minimzation every 1 step is too extreme, I would try with 150 or 250 steps and use the fix box relax along with cg minimization with the vmax option set to 0.0001 to avoid sudden changes in the box volume and convergence problems.
3. The timestep is very short actually. That would not be related to your problem but it would save computation time to increase the timestep.
4. I'm pretty sure that the error you get regarding the memory is due to the very small tdamp in the nvt command you can try the value of 100.0. Because in real units the time is measured in fs and the damping parameter you chose would probably cause the simulation box to explode or yield unphysical temperature values. Probably you can see that if you set your thermo to 1, that if any run was successful with such small damping parameter. In fact, you can try lower than 100, it depends on the system, but I suggest you start from a large value and start to decrease until you get minimum fluctuations in temperature and a stable system. 0.01 would be suitable in metal units not real units.
5. The keyword in the pair_style should be mincap, not minicap. A small typo there.
Thanks again for the detailed explanation. I tried with the minimization with box/relax fix and the volume of the cell kept shrinking to the point i had bonds error, although i did not set any value for vmax. I will try it and see if it changes anything. I will also try with little larger timestep. for the tdamp and pdamp, i am using 100dt but since my dt is very small then those two become small too, once i increase my timestep, i will see what kind of change will have. My temperature was very stable , however my pressure was fluctuating alot even with this small pdamp. miniacap !!! I was certain that i copied it from the website, i guess not :D no wonder lammps wouldnt accept it although lamps log file very low in details.
Thanks again for your help, hopefully i will get it running soon.
- Badges
- Science topic
- Similar topics
- Computer Communications (Networks)
- Simulators