This is an error that occurs when the number of neighbors on any atom exceeds 24 when determining molecular species (in fix reax/c/species). When this error occurs it usually means you have a system that is rapidly changing. From your description, I believe your system is significantly changing with fix temp/berendsen. Visualize the system and the thermodynamic variables to see what has happened.
There is the MAXSPECBOND in reaxc_defs.h but changing it won’t help your simulation too much.
Ray
Hi , thanx for your answer . I was running simulations and playing a bit, and I saw that during nve+temp/berendesen temperature was jumping very high from 300K to 3000K etc and I guess thats why I as getting that mistake of : Increase MAXSPECBOND in reaxc.defs.h .
After that I was just running nve , and after a while NVT or nve+temp/berendesen and I was not getting that mistake. I also moticed that NVT was controlling temperature better then nve+temp/berendesen and speed was quite similar .
However , after running system for while in NVE and NVT I have tried to run it in NPT and I have got again same mistake : : Increase MAXSPECBOND in reaxc.defs.h . After that I have used read_restart and I tried to run simulation in NVE and NVT for longer time and then I have switched to NPT , but again I have got same mistake after I have tried with NPT . I dont understand why is this mistake appearing and what should I do to run simulation in NPT , try shorter timestep or what ? (now I use timestep 0.1 and density is arround 0,3 while equilibrated debsity should be arround 0,8) .
Also what I dont understand about fix reax/c species , when I read species file , it is showing me many small fragments and radicals , and not the molecules which I should have in simulation, does it mean that I should average spieces after every 1000 or 10 000 timesteps (instead if every 100) or after equilibration system will “form” molecules which should be there ?
Also what I dont understand about reax is that I need to spend some time in equilibrating system , but reactions can happen during equilibration too , does it mean that I should equilibrate system on lower temperature and then simulate reaction on higher tempeature , or in reax reactions are reversible so after equilibration on same temperature I will get “average” reaction rate ?
I have one more issue with simulation , my solvent is DMSO , but double bond C=S is not parametrized in reax force fields (but people from reax centre have told me to replace it with some another C=S bond from another force field) , but when I try to run simulation I imediately get either mistake : : Increase MAXSPECBOND in reaxc.defs.h or imediately in first timestep that pressure is nan and then I get message “Program.exe has stopped working” and simulation crashes.
If by pair_style reax command I replace S with C I get “acetone” (like alternative solvent) and simulation starts running until NPT (which is the case I already described above) , so any advices why this happens with DMSO ? I guess there is some problem with C=S bond maybe parameters are not good and forces are too big and it says that pressure is nan ?
Any info or explanations or advices would be welcome as I am trying to learn how to run simulations with reax .
Thank you.
Have you visualized your system during/after NVE+temp/berendsen and NPT runs? This error is associated with atoms having more than 24 bonds – which is an indication of bad/weird dynamics.
How are you using the fix reax/c/species command? Did you use time-averaging or the “cutoff” keyword? From the doc page, “Optional keyword cutoff can be assigned to change the minimum bond-order values used in identifying chemical bonds between pairs of atoms. Bond-order cutoffs should be carefully chosen, as bond-order cutoffs that are too small may include too many bonds (which will result in an error), while cutoffs that are too large will result in fragmented molecules. The default cutoff of 0.3 usually gives good results.”. But apparently the default value of 0.3 didn’t work for you. You would need to experiment with the cutoffs.
Reactions can take place at any time and at any temperature if it more energetically favorable. There are no definite reversabilities, broken bonds can re-form and existing/new bonds can break. You can monitor chemical reactions with fix reax/c/species right from the beginning of the equilibration.
Are you sure the particular ReaxFF description you are using is capable of describing C and S with a double bond? What the “people from reax center” said does not make sense to me. You cannot mix-and-match ReaxFF descriptions. If the simulation fails with S but works with C, that just means the parameter set you used is not suitable for S species.
Note that there is no C-S, C=S, or C triply-bonded to S, per se. ReaxFF determines bond order, strength, and energy on-the-fly. You don’t specify it to be “double” and expect the bond to stay that way.
Ray
I also think this error is misleading. Simply increasing MAXSPECBOND won’t work since the number 24 is hardcoded into compute_spec_atom.cpp.
Stan
Would it make more sense to just terminate the run (with an explanatory error code) if this condition occurs?
Jim
Hi , thanx for the explanations, I still didnt visualize system yet but I will do it during next runs and try to increase bond cutoff and see how will simulation go (and weird thing is that when I run same simulation on server instead on computer I dont get this error) .
I have wrote to Reax research center whose link was given on LAMMPS website and they told me to
potentially, transfer the S-O interaction parameters from
Shin, Y.K., Kwak, H., Vasenkov, A., Sengupta, D. and van Duin, A.C.T. (2015) Development of a ReaxFF reactive force field for Fe/Cr/O/S and application to oxidation of butane over a pyrite-covered Cr2O3 catalyst. ACS Catalysis 5, 7226-7236.
to the protein force field - and this may work for a non-reactive DMSO.
So basically if simulation works with C and not with S it is because of bad parameters and it has nothing to do with LAMMPs , I should find another “alternative” solvent or maybe ask people from reax about it .
About reax/c species I still dont understand that I am getting many fragmented molecules or some very ilogical molecules (for example CaromCdmso14CCHd1HmK2Odmso3Cs6) but maybe increasing bond cutoff and running simulation for longer time could help.
Hi , I am trying to use cutoff of keyword and I am getting error : Ilegal fix ave/atom command
(for command fix a all reax/c/species 1 1000 1000 species.out cutoff 2 3 0.5 element Ca C C C Ha H H H H K N O O C ) , if I use fix a all reax/c/species 1 1000 1000 species.out cutoff * * 0.5 element Ca C C C Ha H H H H K N O O C
I am getting error :
illegal fix reax/c species command .
and I am using same command like from LAMMPS example page , I dont understand what is happening
Fix reax/c/species’ cutoff keyword does not allow wildcard (I.e., asterisk or *) in its argument – as it was not advertised on the doc page.
Ray
Yes but I said that I have also used command :
fix a all reax/c/species 1 1000 1000 species.out cutoff 2 3 0.5 element Ca C C C Ha H H H H K N O O C
and its still giving mistake Ilegal fix ave/atom command .
Worked fine for me. Could you try with the latest lammps version and if the error persits attach the full input deck? At this stage it is impossible to move on without your input files.
Ray
Hi , I am using latest version of LAMMPS ICMS 27th August, I have attached the log file from simulation and error it gives me .
I wanted to try simulation with larger cutoff but I am getting this mistake .
log.lammps (1.02 KB)
I can’t reproduce your simulation with a log file…
The only thing I can do is to suggest you re-write the fix reax/c/species command as the failure came from “1 1000 1000” - there might be texting issues.
Ray