[gmx-users] equilibrating a solvent

[Sidney Elmer]

Erik Lindahl wrote:

> Sidney Elmer wrote:
>
> >I would appreciate any help anyone could give me on equilibrating a
> >carbon tetrachloride solvent I made.  I've created the topology file and
> >the .gro file, now I want to equilibrate it in an NVT ensemble.  What
> >I've done is written a script that will create a
> >crystal of the solvent.  I have checked to make sure that I've placed
> >the necessary number of molecules in the box of a given volume to give
> >the desired experimental density.  Now I just have to equilibrate the
> >liquid.  In the .mdp file, whether I specify to use md with Berendson
> >temperature coupling or stochastic dynamics, I get the same error
> >message from grompp, which
> >I have attached to this email.  There are 2560 atoms in the system, but
> >for some reason, it appears that not all of the atoms are coupled to the
> >
> >heat bath.  As a side note, I noticed the first time I ran grompp, the
> >error message specified that it was "Making dummy/rest group for
> >T-Coupling containing 2555 elements" instead of the 2560.  That is why
> >I'm thinking that some of the atoms aren't being coupled to the
> >temperature bath, (but it may be unrelated).  I just don't understand
> >the Fatal error:  Not enough ref_t and tau_t values.  I will also
> >include the .mdp file for you as well.  I won't bother sending the .gro
> >file, because it's pretty big, and it's sufficient to know that there
> >are 2560 atoms in the file, and the box size is correct on the last
> >line.  If you could help me I'd appreciate it.  Thanks.
> >
> >Sid
> >
> Hi Sid,
>
> The first thing to do is probably an energy minimization. Use
> integrator=steep
> and run say 100 steps, just to get rid on any atom overlaps. Turn off
> the constraints
> during the energy minimization.
>
> Once you have a reasonable (i.e., negative) energy, you can run SD or MD
> with
> temperature coupling. The temperature coupling stuff seems to be
> commented out in your mdp file.
> Try using e.g:
>
> tcoupl = Berendsen
> tc-grps = system
> tau_t = 0.1
> ref_t = 300
>
> And turn of the pressure coupling to start with (keeping it commented
> out will do that). Once
> you have equilibrated the box for a little while you can turn on
> pressure coupling too...
>
> Cheers,
>
> Erik

Hi Erik,

I did the energy minimization, that worked fine.  But the problem is with
grompp.  I did comment out the Berendson temperature coupling in the .mdp
file because I wanted to use Stochastic Dynamics, so I set integrator to
"sd".  Even though the integrator is set to "sd", I still get the
Error Message that there are not enough ref_t and tau_t values.  So I assume
there is a similar mechanism for coupling atoms to a temperature bath,
whether I use stochastic dynamics or Berendson temperature coupling.  Either
way, I get the same error.

I only have one group definition, CCL4, so making a list of values in .mdp
file for ref_t and tau_t didn't work.  I am confused by the error, not so
much by the sequence of events needed to get an equilibrated system.  I also
have my topology file, maybe there is a problem there?
-------------- next part --------------
;     TESTING CCL4 TOPOLOGY FILE
;     WRITTEN BY SID ELMER
;     JUNE 13, 2002
;
;  Include forcefield parameters
#include "ffG43a1.itp"

[ moleculetype ]
; name  nrexcl
CCL4	2

[ atoms ]
;   nr    type   resnr  residu  atom    cgnr    charge   mass
    1     CCL4   1	CTC 	CCl 	1	 0.248   12.011
    2     CLCL4  1	CTC 	CL1  	1	-0.062   35.453
    3     CLCL4  1	CTC 	CL2  	1	-0.062   35.453
    4     CLCL4  1	CTC 	CL3  	1	-0.062   35.453
    5     CLCL4  1	CTC 	CL4  	1	-0.062   35.453

[ bonds ]
;  ai    aj funct          b0       kb
    1     2     1 	0.1769	8.1000e+06	; parameters taken from Jorgensen
    1     3     1 	0.1769	8.1000e+06	; JACS, 114 (19), pp7535, 1992
    1     4     1 	0.1769	8.1000e+06	; reference #18
    1     5     1 	0.1769	8.1000e+06	; 
;    2     3     1 	0.2889	3.0400e+06	; subsequently converted to gmx units
;    2     4     1 	0.2889	3.0400e+06
;    2     5     1 	0.2889	3.0400e+06	; only r(c-cl) given, r(cl-cl) 
;    3     4     1 	0.2889	3.0400e+06	; calculated by geometry.
;    3     5     1 	0.2889	3.0400e+06	; but no kb given for either bond type,
;;    4     5     1 	0.2889	3.0400e+06	; due to rigid model.  Therefore, guess.

[ angles ]
;    i     j      k     funct      a0        ka
     2     1      3       1       109.5    618.00
     2     1      4       1       109.5    618.00
     2     1      5       1       109.5    618.00
     3     1      4       1       109.5    618.00
     3     1      5       1       109.5    618.00
     4     1      5       1       109.5    618.00

;[ constraints ]
;;   i     j    funct      dij
;    1     2      1       0.1769
;    1     3      1       0.1769
;    1     4      1       0.1769
;    1     5      1       0.1769
;    2     3      1       0.2889
;    2     4      1       0.2889
;    2     5      1       0.2889
;    3     4      1       0.2889
;    3     5      1       0.2889

[ exclusions ]
;   i     j     k
    2     1     3
    2     1     4
    2     1     5
    3     1     4
    3     1     5
    4     1     5


[ system ]
; Name
CCL4 single molecule

[ molecules ]
; Compound      #mols
CCL4             512