When uniaxially compress loads of cnt (along x direction, under nvt, cnt could buckl), I use fe=pxx*ly*lz to calculate the external force.
The problem is that some portion of the atoms around x=0 are frozen through fix setforce 0 0 0, which leads to the drop of obtained fe, as posted in this thread. Of course, I’ll use fix rigid instead in the future. Nevertheless, is it possible to re-gain the realistic pressure after using setforce 0 0 0 ? The simulation is quite time-consuming.
I have read the compute pressure in the manual, and understand that pxx is calculated as (only the component pxx is needed)
pxx=1/V*sum(m_k*v_k^2,{k,1,N})+1/V*sum(r_k*f_k,{k,1,N’})
In my simulation, T=0.1K and the nominal strain ratio (of the box)=1e-5/ps, so the first term may be not so important. Two additional questions rise: a) how to determine r_k? if r_k is the x coordinate, then it seems that shifting the box would change the results. b) how to determine N’ ?
Not really. Bascially, you have an inconsistent system and there is no clean way out of it.
You could use compute stress/atom on the atoms that are let mobile and compute the sum over those, but then you need know the exact volume to convert the virial stress to pressure and that volume is not well defined but must be approximated. This gives you the pressure only for the volume of the mobile atoms and it includes the contributions from the forces of the immobile atoms to the mobile atoms but not vice versa. So this is not 100% consistent and clean.
You could use neigh_modify exclude to exclude the interactions within the the immobile block of atoms, but have them otherwise contribute to the global pressure. But then the kinetic contribution of the pressure also includes the immobile atoms and setting that to 0 is not 100% consistent and clean.
You can use fix rigid on the immobile atoms and have a correction for the pressure due to keeping a rigid body, but this also does not account for the “temperature” of the rigid body, which is set to an unphysical 0K (or worse some arbitrary value).
Instead of using fix setforce 0 0 0, you could use fix spring/self. This way the immobile atoms are not fully immobile but can have interactions and move as far as the attached spring allows. That avoids the inconsistencies from above but now includes an artificial restraining force. How well, this works also depends on how well relaxed that immobile group of atoms is and thus how much artificial force is needed to sufficiently restrain those atoms.
If you have a trajectory with positions and velocities, you can experiment with the different settings on the same trajectory data using the rerun command.
Bottom line, however you turn it, you cannot avoid some degree of inconsistency; the best you can do is figure out what minimizes the impact for your property of interest.
Thank you so much @akohlmey for your systematic and constructive suggestions and analyses. It provides valuable guidance to my coming simulations.
I guess it would not be possible to re-gain the reasonable approximation of fe=pxx*ly*lz through the already done simulations, which used fix setforce.
I have another question. Suppose both the temperature and loading speed are sufficiently low, so the kinetic could be neglected. Is the effect of ‘setforce 0 0 0’ on fe constant? I mean is it possible that there is a certain unkonwn constant k during the whole simulation, such that k*pxx*ly*lz becomes a reasonable approximation of the external compression loads when ‘setforce 0 0 0’ is used? Please note that cnts may buckle during compression, so the deformation along x direction may not be uniform.
My previous post was based on general MD knowledge and observations. You are now asking about something that is specific to your research and outside my area of expertise in general MD. So I cannot provide a meaningful answer.