[gmx-developers] pbc and virial problem

Wed Apr 8 17:09:38 CEST 2009

Hi,

I still can't make much of this.

Is is not simply a incorrect memory access error, where you write to 
fshift iso to the forces?
You can run valgrind to check for that.

Alexander Herz wrote:
> 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)).
>   
No, these are shifts by 1 box vector, not 3.

Berk
> 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.
>>   
>>     
>
> _______________________________________________
> 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.
>   