[gmx-users] ewald_geometry 3dc with wall gives high pressure

[zifeng li]

>> Zifeng Li lizifeng7 at gmail.com
>> Dear Gromacs users,
>>
>> I want to compress a polymer slab and solvate it with water to exam the interface between them. To prevent the polymer from crossing the interface along z direction, I compress it using 2 walls by moving them towards each other a little bit ( 0.5 sigma of wall atoms) and equilibrating for 50ps subsequently each time. However, the pressure is enormously high during this process.
>>
>> From polymer handbook, this slab should have a density of 0.98g/cm3 at 400K, 1bar.
>> In my simulation, the density-pressure relation is :
>>
>> density(g/cm3)     Avg. pressure(/bar)
>> 0.45                     80
>> 0.63                     245
>> 0.86                     2500
>>
>> Gromacs version: 4.5.4
>>
>> Here's my NVT.mdp settings:
>>
>> constraints         =  none
>> integrator          =  md
>> dt                  =  0.001; 1fs.
>> nsteps              =  50000; total 50ps.
>> ;neighbour searching
>> nstlist             =  10
>> rlist               =  1.1
>> ns_type             =  grid
>>
>> ;coulomb
>> coulombtype       =  PME
>> rcoulomb            =  1.1
>> ewald_geometry  = 3dc
>>
>> ;walls(9-3 type by default)
>> pbc                = xy
>> nwall =2
>> wall_atomtype          =opls_139 opls_139
>> wall_density          = 50 50  (  set is as 1g/cm3)
>>
>> ;vdw
>> rvdw                = 1.1
>> vdwtype             = shift
>>
>> ;output control
>> nstxout             =  1000
>> nstvout             =  0
>> nstfout             =  0
>> nstlog              =  1000
>> nstenergy           = 1000
>> ; Generate velocites is on at 400 K.
>> gen_vel             =  yes
>> gen_temp            =  400.0
>> gen_seed            =  173529
>>
>> ;temperature coupling
>> tcoupl                   = v-rescale
>> tc-grps                  = Polymer
>> tau_t                    = 0.1
>> ref_t                    = 400
>>
>> equilibrating the output of compression which has density 1.0 and run a NPT with walls, it keeps expanding to density 0.1 g/cm3 after 6ns.
>> I doubt this high pressure issues might be caused by wall settings or ewald_geometry settings.
>>
>> 1. For the wall settings
>>
>>  The density of wall (1g/cm3) is close to the density of my system (0.45~1g/cm3) and each time I run a NVT, the distance between wall and system is 0.5 sigma of wall atoms, the energy between wall and system is 4% of the total energy of my system.  I don't know what's wrong or what I can check next.

       any problem with wall settings?
>
>
>>Therefore, I just remove the wall and check the ewald_geometry.

>> 2. For ewald_geometry = 3dc
>>
>> >From the paper J.Chem.Phys.111(1999)  pp3155-3162 on Ewald summation for systems with slab geometry, it requires the length of box in z direction should be at least 3 times of that in x/y direction to use this psudeo-2D PME.
>>
>> However, to simulate such a thick slab would be computationally expensive. Therefore, I put vacuum above and below the slab(at slab density 0.45), running a NVT without walls(10ns), the pressure goes down to some negative value near 0 and the polymer compresses itself in z direction to density avg. 0.63 along this NVT and stops shrinking. The density(0.63) is too low compared to my expectation(1.0g/cm3). However, I check the density profile and find in the middle of the slab, it has bulk density(nearly 1.0g/cm3) and the density decrease towards the interface. It seems reasonable since it has right bulk density in the middle but I lose control of pressure( along z, it is 0 because of vacuum) and density( the avg. density value is low).
>
>
>    Is this a correct or efficient way to use 3d ewald correction when you don't want a thick slab?
>>
>>
>> Thanks!
>>
>> --Zifeng