[gmx-developers] pbc and virial problem

[Alexander Herz]

Ok,

I still don't understand why this happens, but I've got some more
details now.
I have averaged the different virial contributions over 50*10^3 steps:

With correction:
50000:
AV_HOME_VIR:[-39226.070555,-39294.997520,-37501.780167]
//just [f[i][XX]*x[i][XX],f[i][YY]*x[i][YY],f[i][ZZ]*x[i][ZZ]] over all
home atoms

AV_SHIFT_VIR: [594.654476,582.680015,0.000641]
//same for fshift only


fshift[7]={-0.000129,-0.000005};fshift[8]={-0.000079,0.000003};fshift[16]={42.458567,38.451720};fshift[17]={64.619434,56.483254};fshift[18]={-38.068221,38.650424};fshift[21]={226.567056,-32.808623};fshift[22]={-1510.309930,-302.387872};fshift[23]={-28.218844,-56.201202};fshift[26]={20.397865,-50.388570};fshift[27]={96.795412,-180.068616};fshift[28]={-40.725687,-24.410759};
//fshift[0..SHIFTS] averaged (x and y only)
av_fshift={-1189.308714,-1165.360030};
//AV_SHIFT_VIR*-2 (calculated from the averaged fshifts above)

AV_EL_RECIP_VIR:[-4621.878257,-4635.251526,-4433.260378]
AV_SHAKE_VIR:[47724.981107,47806.627009,46467.525646]
AV_TOTAL_VIR:[4471.686771,4459.057978,4532.550496]

Without correction:
50000:
AV_HOME_VIR:[-39201.051820,-39234.433725,-37260.849326]
AV_SHIFT_VIR:[942.510704,943.765406,0.000000]
fshift[16]={40.338618,46.695308};fshift[17]={-23.836940,215.274137};fshift[18]={-42.609950,50.021658};fshift[21]={165.046904,32.068133};fshift[22]={-118.993265,-668.705978};fshift[23]={-293.527576,16.681827};fshift[26]={53.153976,-33.448852};fshift[27]={-23.787922,-259.604220};fshift[28]={-33.663446,-24.132762};
av_fshift={-1885.021409,-1887.530811};
AV_EL_RECIP_VIR:[-4621.740191,-4630.363645,-4398.286901]
AV_SHAKE_VIR:[47690.282258,47747.270976,46197.490017]
AV_TOTAL_VIR:[4810.000951,4826.239012,4538.353789]

The main difference in the total virial clearly can be attributed to the
shift virial (AV_SHIFT_VIR calculated from fshift[0..SHIFTS]). Looking
at the different averaged fshift contributions, fshift[17,21,23,27]
contribute most
(these correspond to the shifts (0,-3,0),(-3,0,0),(3,0,0),(0,3,0)).

I'm using spc/e water by the way, the correction force will act on the
OW atoms only.

Any ideas why the shifted virial contributions differ so much while the
unshifted (home virial) contributions
are pretty much unaffected (as one would assume for all contributions in
x/y)??

How are the shifted forces calculated? Is there a specific
optimization/assumption for the
spc/e water optimized kernels I might be breaking? Aren't the shifted
forces depending on the "home" forces?
I really don't get why the shifted once should differ when the home ones
don't??

I'm really a bit lost here.

Thx,
Alex


hessb at mpip-mainz.mpg.de schrieb:
> Hi,
>
> What you are saying seems a bit contradictory.
> Do I understand it correctly that you only add something
> to the z component of the force (and not to the shift force)
> and then the components of the shift forces change, which
> affect only the xx and yy components of the virial?
> That seems very strange.
> This almost must mean that somehow the position of atoms
> along x and y changes due to your extra force.
>
> For your force you have taken care of the correction that was missing
> in the wall potentials that caused the asymmetry?
>
> Berk
>
>   
>> Hi,
>>
>> I've added a force similar to the walls. Now I'm trying to calculate
>> surface tensions using this new correction via the virial.
>> We used to do these calculations with pbc=xyz with an older gromacs
>> version so I'm trying to do the same in order to be able compare the
>> results before going to pbc=xy.
>>
>> The problem is that I get non sense for the virial xx and yy components.
>> To figure out why, I compared the averaged virials for the same sim with
>> and without the new correction.
>>
>> The virial contribution from the non-shifted particles is pretty much
>> identical but the virial from the shifted
>> particles:
>>   (calc_virial in sim_util.c):
>>   clear_mat(vir_part);
>>
>> calc_vir(fplog,SHIFTS,fr->shift_vec,fr->fshift,vir_part,ePBC==epbcSCREW,box);
>>   inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
>>
>> differs by almost exactly the amount which I'd expect.
>>
>> The correction I apply only adds forces in z direction and the
>> simulation is setup such that no pbc in z is applied because of the cut
>> off (there is a water slab of ca. 6 nm thickness in the middle of the
>> box (z direction) surrounded by ca. 3 nm of nothing in both z
>> directions). Debugging through the code I can see that no images for
>> z!=0 are created.
>> The correction force is quite small compared to the forces applied by
>> gromacs non_bonded
>> (correction: avrg. force=0.03, max. force=8.0, gromacs: avrg.:456.0 ,
>> max.: 4831.0).
>>
>> The force is calculated just after the walls right before do_nonbonded.
>>
>> Finally, the question:
>> Any ideas how the small additional force in z direction can have such a
>> strong effect on the virial xx and yy components
>> (apparently fshifted[XX] and fshifted[YY] are influenced significantly
>> whereas the virial ZZ value remains virtually identical
>> no matter if the correction is applied or not, as fshifted[ZZ]==0 in
>> both cases).
>>
>> Thx,
>> Alex
>>
>> _______________________________________________
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://www.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-developers-request at gromacs.org.
>>
>>     
>
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.
>