[gmx-users] Implicit solvent models in Gromacs ?

[Graham Smith]

Hi Marc, 

I think the explanation you want is on p 113-114 of the manual. 
(section 6.4). It was also discussed last month on gmx-users I think?

A bit of code  (which hasn't been properly
tested!) to do the epsilon = r-in-Angstoms thing is 




# include <stdio.h> 
# include <math.h> 
# include <string.h> 

# ifdef DOUBLE 
# define DR 0.0005 
# else 
# define DR 0.002 
# endif 

void get_power(float r, float *fr, float *d2fr, int n) 
{
  if (r<0.039999) {
    *fr = 0.0 ; 
    *d2fr = 0.0 ; 
  }else{
    *fr = exp(n*log(r)) ; 
    *d2fr = n*(n-1)*exp((n-2)*log(r)) ; 
  }
}


int main () 
{
  int i, j ; 
  float max = 2.00001 ; 
  int nstep = max/DR ; 
  float r ; 
  float f1, d2f1 ; 
  float f2, d2f2 ; 
  float f3, d2f3 ; 
  FILE *OUT ;  
  char fname[]="table-2-6-12.xvg" ; 
  
  OUT=fopen(fname,"w") ; 
  for (i=0 ; i<=nstep ; i++ ) {
    r = i * DR ; 
    get_power(r, &f1, &d2f1, -2) ; 
    f1 *= 0.1 ; d2f1 *= 0.1 ; 
    get_power(r, &f2, &d2f2, -6) ;  
    f2 *= -1 ; d2f2 *= -1 ; 
    get_power(r, &f3, &d2f3, -12) ; 
    fprintf(OUT, "%10.3f  %10.6g  %10.6g  %10.6g  %10.6g  %10.6g
%10.6g\n", 
	    r,   f1, d2f1, f2, d2f2, f3, d2f3 ) ; 
  }
}


########################################################################

Dr. Graham R. Smith, 
Biomolecular Modelling Laboratory, 
Imperial Cancer Research Fund, 
44 Lincoln's Inn Fields, 
London WC2A 3PX,
U.K. 
Tel: +44-(0)20 7269 3348 
email: graham.smith at icrf.icnet.uk
URL: http://www.bmm.icnet.uk/~smithgr

On Wed, 12 Dec 2001, Marc Baaden wrote:

> 
> Sorry, .. I first sent the reply to David only, not to the list ...
> 
> >>> David van der Spoel said:
>  >> You have to make a user defined potential (using tables).
>  >> You can then use V=1/r^2 for the Coulomb if you so wish, but check your
>  >> units! People usually have this when r is in Angrstrom, so in the GROMACS
>  >> case you need a factor 10 somewhere. Check the file table6-12.xvg in the
>  >> top directory, that you have to modify, and then pass it to mdrun.
> 
> OK. That is fine.
> BUT .. what format is that table file ? From the handbook page 166 I
> would have expected 4 columns yi, Fi, Gi and Hi, but there are 7 of those ..
> 
> 
> Also I wonder whether there is any example program/script, ..
> which shows how table6-12.xvg is calculated (and even better
> where I could just replace the epsilon-r by 1/r^2 to get what
> I want ....)
> 
> 
> Which is the actual part (function) stored in the table ?
> E.g. if you have
> 
>                qi * qj
> Vcoulomb = f * --------
>                er * rij
> 
> I suppose the tabulated function would be
>    
>            f
> y(x) = --------
>         er * x
> 
> where x=rij
> 
> as qi and qj depend on the atom types of at i and j and cannot be included
> in the table. Is that right ?
> 
> Am I also right that you actually have to store values only up to the
> Cutoff you are going to use (eg the values in table6-12 go from 0 to
> 2 (nm), which should be more than enough).
> 
> 
> What are the files table6-10.xvg, table6-11.xvg, table6-8.xvg and
> table6-9.xvg ?
> 
> 
> I also wonder wether there are any example ctab.xvg, rtab.xvg, .. files
> for Gromacs 2.0. They seem to be of a different format (Gromacs complains
> and says "Fatal error: Trying to read file ctab.xvg, but no colums = 7,
> should be 5 "
> 
> 
> Marc
> 
> -- 
>  Dr. Marc Baaden - Laboratory of Molecular Biophysics, Oxford University
>  mailto:baaden at smplinux.de  - ICQ# 11466242 -  http://www.marc-baaden.de
>  FAX/Voice +49 697912 39550  -  Tel: +44 1865 275380  or  +33 609 843217
> 
> 
> 
> -- 
>  Dr. Marc Baaden - Laboratory of Molecular Biophysics, Oxford University
>  mailto:baaden at smplinux.de  - ICQ# 11466242 -  http://www.marc-baaden.de
>  FAX/Voice +49 697912 39550  -  Tel: +44 1865 275380  or  +33 609 843217
> 
> 
> 
> -- 
>  Dr. Marc Baaden - Laboratory of Molecular Biophysics, Oxford University
>  mailto:baaden at smplinux.de  - ICQ# 11466242 -  http://www.marc-baaden.de
>  FAX/Voice +49 697912 39550  -  Tel: +44 1865 275380  or  +33 609 843217
> 
> 
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.
>