On Wed, Mar 9, 2011 at 8:54 AM, David A Case <case.biomaps.rutgers.edu> wrote:
> On Mon, Mar 07, 2011, Andrew Jewett wrote:
>
>> If the amber community is considering
>> adding new features to AMBER, then please consider the following:
>>
>> i) The ability for a covalent bond to cross the periodic box boundary.
>
> This is likely to happen.
Hello.
I'm glad to hear that. That will be useful.
>> ii) The ability to allow pressure rescaling in only the Z direction.
>> (while preventing rescaling in the XY directions, or visa versa).
>
> This sounds almost like a two-line change to the code -- have you considered
> trying it out?
Yes.
(For the record, my tests, and also the figure with the bubble use
this. But it's hard-coded.)
-------------------------
In case somebody else wants to try this out, here is the code change
you need to make (at least in sander). (If it's okay, I can post the
whole files.) Here are the changes. Go to "runmd.f" and replace the
following lines of code (near line 2390) :
else if (ntp > 1) then
rmu(1) = (1.d0-dtcp*(pres0-ener%pres(1)))**third
rmu(2) = (1.d0-dtcp*(pres0-ener%pres(2)))**third
rmu(3) = (1.d0-dtcp*(pres0-ener%pres(3)))**third
end if
---replace with:---
else if (ntp > 1) then
rmu(1) = 1.d0
rmu(2) = 1.d0
rmu(3) = (1.d0-dtcp*(pres0-ener%pres(3)))**third
end if
>> a) Harmonic restraint anchor points are not fixed during NPT simulations,
>> but are scaled. It would be useful to make this optional
>> (both for belly and harmonic restraints).
>
> See comment above.
Sure. It's not hard, and I'll post the code below for those interested.
But I strongly suspect that changing this code alone won't prevent the
bubble problem. (...But I'll test this and get back to you today.) If
I understand things correctly, the bubble problem is a consequence of
the way the virial is calculated, (and the fact that I turned off the
charge). I think this change needs to be accompanied by a
modification to the virial.
-------------------------
For the record, here's the changes I made to the code. (Disclaimer:
I'm very new to fortran.)
By the way, I'm using belly to constrain the atom positions. (I like belly.)
(If you want to stick with harmonic restraints, alas,
you probably have to make modifications to the code somewhere else.)
Go to ew_box.f and make three changes to the code:
1) near line 27 of ew_box.f, I added an argument, irescalemask, which
is an array of integers to subroutine ew_pscale()
Here's the corresponding code change. Replace:
subroutine ew_pscale(natom,x,amass,nummols,molsiz,npscal)
use nblist, only: oldrecip,ucell
implicit none
---with:---
subroutine ew_pscale(natom,x,amass,nummols,molsiz,npscal, &
ibelly,irescalemask)
use nblist, only: oldrecip,ucell
implicit none
integer ibelly,irescalemask(*) ! <--I added this line
2) Near line 89 of ew_box.f, replace the following do loop:
! ...now rigidly translate molecule
do iatom = 1,num
i = i + 1
x(1,i) = x(1,i) + xmolnu - xmol
x(2,i) = x(2,i) + ymolnu - ymol
x(3,i) = x(3,i) + zmolnu - zmol
end do
---with:---
! ...now rigidly translate molecule
do iatom = 1,num
i = i + 1
if ((ibelly /= 2).or.(irescalemask(n) /= 0)) then
x(1,i) = x(1,i) + xmolnu - xmol
x(2,i) = x(2,i) + ymolnu - ymol
x(3,i) = x(3,i) + zmolnu - zmol
end if
end do
This way, whether or not the belly atoms are rescaled remains optional.
When ibelly=1, the belly atoms are rescaled, which is the original behavior.
When ibelly=2, the belly atoms are not rescaled.
3) Make a similar change in ew_box.f near line 118. Replace:
! ...use these with the new cell parameters to get new
! cartesian coordinates
x(1,n) = f1*ucell(1,1)+f2*ucell(1,2)+f3*ucell(1,3)
x(2,n) = f1*ucell(2,1)+f2*ucell(2,2)+f3*ucell(2,3)
x(3,n) = f1*ucell(3,1)+f2*ucell(3,2)+f3*ucell(3,3)
---with:---
! ...use these with the new cell parameters to get new
! cartesian coordinates
if ((ibelly /= 2).or.(irescalemask(n) /= 0)) then
x(1,n) = f1*ucell(1,1)+f2*ucell(1,2)+f3*ucell(1,3)
x(2,n) = f1*ucell(2,1)+f2*ucell(2,2)+f3*ucell(2,3)
x(3,n) = f1*ucell(3,1)+f2*ucell(3,2)+f3*ucell(3,3)
end if
4) Finally, go to "runmd.f", and modify every place you invoke
ew_pscale() accordingly. For example, replace
call ew_pscale(natom,x,amass,nspm,nsp,2)
---with:---
call ew_pscale(natom,x,amass,nspm,nsp,2,ibelly,ix(ibellygp))
and later, replace:
call ew_pscale(natom,x,amass,nspm,nsp,npscal)
---with:---
call ew_pscale(natom,x,amass,nspm,nsp,npscal,ibelly,ix(ibellygp))
>> b) Arguably, constraints cause the virial to be calculated incorrectly.
>> This can force the box to expand to a huge size.
>
> What kind of "constraints" are you using? Do you understand what is going on,
> e.g. why there is a bubble only on one side of the central slab?
Yeah, I was surprised too when I saw that, but it happened over and
over again in different test runs. I should note, we weren't using a
very realistic ice structure. In our ice crystal, all the hydrogen
bond acceptors are on the top basal plane of the ice crystal, and the
donors, are on the bottom face. That's probably related to why the
bubble always appeared on the top face. It seems odd to me too.
Incidentally, in my last post I forgot to include a caption for the
bubble figure.
Here are the details for the bubble figure (figure caption):
"Here we show a frame animation from a simulation of (default TIP3P)
water interacting with a slab of immobilised ice (middle), 6 water
molecules thick, which has been held in place using harmonic
restraints under NPT conditions (T=300K,$P=1 bar). The inner 4 layers
of ice (white) have had their charges removed (by editing the .top
file), allowing the remaining Lennard-Jones forces to dominate.
However, the outermost layer of water on either side of the ice slab
is unmodified."
SIDE-COMMENT: I think the Lennard-Jones radii of the atoms in TIP3
water molecules is actually fairly large... larger than the spacing
between them. Once the charges are removed, we think there is
significant intermolecular repulsion between the atoms.
"Under NPT conditions, AMBER continues to expand the simulation box in
a vain effort to relax the stress between repulsive (white) molecules,
which are held in place by restraints. This cavitation problem goes
away if either:
1) the charges in the middle layer are restored, or
2) harmonic restraints are removed allowing the neutral (white)
molecules to expand to their equilibrium density."
----As for why somebody would run a simulation like this:---
A number of simulations of neutralised, immobilised water have been
carried out by the groups of Shekhar Garde and David Chandler and
others to study the hydrophobic effect, although I think they use SPC
water.
Cheers!
Andrew
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Mar 09 2011 - 13:30:04 PST
