src/tools/gmx_membed.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2012, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <signal.h>
  43 #include <stdlib.h>
  44 #include "typedefs.h"
  45 #include "sysstuff.h"
  46 #include "statutil.h"
  47 #include "macros.h"
  48 #include "copyrite.h"
  49 #include "main.h"
  50 #include "gmx_ana.h"
  51 #include "string2.h"
  52
  53 int gmx_membed(int argc,char *argv[])
  54 {
  55     const char *desc[] = {
  56         "[TT]g_membed[tt] embeds a membrane protein into an equilibrated lipid bilayer at the position",
  57         "and orientation specified by the user.[PAR]",
  58         "SHORT MANUAL[BR]------------[BR]",
  59         "The user should merge the structure files of the protein and membrane (+solvent), creating a",
  60         "single structure file with the protein overlapping the membrane at the desired position and",
  61         "orientation. The box size is taken from the membrane structure file. The corresponding topology",
  62         "files should also be merged. Consecutively, create a [TT].tpr[tt] file (input for [TT]g_membed[tt]) from these files,"
  63         "with the following options included in the [TT].mdp[tt] file.[BR]",
  64         " - [TT]integrator      = md[tt][BR]",
  65         " - [TT]energygrp       = Protein[tt] (or other group that you want to insert)[BR]",
  66         " - [TT]freezegrps      = Protein[tt][BR]",
  67         " - [TT]freezedim       = Y Y Y[tt][BR]",
  68         " - [TT]energygrp_excl  = Protein Protein[tt][BR]",
  69         "The output is a structure file containing the protein embedded in the membrane. If a topology",
  70         "file is provided, the number of lipid and ",
  71         "solvent molecules will be updated to match the new structure file.[BR]",
  72         "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.[PAR]",
  73         "SHORT METHOD DESCRIPTION[BR]",
  74         "------------------------[BR]",
  75         "1. The protein is resized around its center of mass by a factor [TT]-xy[tt] in the xy-plane",
  76         "(the membrane plane) and a factor [TT]-z[tt] in the [IT]z[it]-direction (if the size of the",
  77         "protein in the z-direction is the same or smaller than the width of the membrane, a",
  78         "[TT]-z[tt] value larger than 1 can prevent that the protein will be enveloped by the lipids).[BR]",
  79         "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
  80         "intraprotein interactions are turned off to prevent numerical issues for small values of [TT]-xy[tt]",
  81         " or [TT]-z[tt][BR]",
  82         "3. One md step is performed.[BR]",
  83         "4. The resize factor ([TT]-xy[tt] or [TT]-z[tt]) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
  84         "protein is resized again around its center of mass. The resize factor for the xy-plane",
  85         "is incremented first. The resize factor for the z-direction is not changed until the [TT]-xy[tt] factor",
  86         "is 1 (thus after [TT]-nxy[tt] iterations).[BR]",
  87         "5. Repeat step 3 and 4 until the protein reaches its original size ([TT]-nxy[tt] + [TT]-nz[tt] iterations).[BR]",
  88         "For a more extensive method description see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.[PAR]",
  89         "NOTE[BR]----[BR]",
  90         " - Protein can be any molecule you want to insert in the membrane.[BR]",
  91         " - It is recommended to perform a short equilibration run after the embedding",
  92         "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174), to re-equilibrate the membrane. Clearly",
  93         "protein equilibration might require longer.[PAR]"
  94     };
  95     t_filenm fnm[] = {
  96         { efTPX, "-f",      "into_mem", ffREAD },
  97         { efNDX, "-n",      "index",    ffOPTRD },
  98         { efTOP, "-p",      "topol",    ffOPTRW },
  99         { efTRN, "-o",      NULL,       ffWRITE },
 100         { efXTC, "-x",      NULL,       ffOPTWR },
 101         { efSTO, "-c",      "membedded",  ffWRITE },
 102         { efEDR, "-e",      "ener",     ffWRITE },
 103         { efDAT, "-dat",    "membed",   ffWRITE }
 104     };
 105 #define NFILE asize(fnm)
 106
 107     /* Command line options ! */
 108     real xy_fac = 0.5;
 109     real xy_max = 1.0;
 110     real z_fac = 1.0;
 111     real z_max = 1.0;
 112     int it_xy = 1000;
 113     int it_z = 0;
 114     real probe_rad = 0.22;
 115     int low_up_rm = 0;
 116     int maxwarn=0;
 117     int pieces=1;
 118     gmx_bool bALLOW_ASYMMETRY=FALSE;
 119     gmx_bool bStart=FALSE;
 120     int nstepout=100;
 121     gmx_bool bVerbose=FALSE;
 122     char *mdrun_path=NULL;
 123
 124 /* arguments relevant to OPENMM only*/
 125 #ifdef GMX_OPENMM
 126     gmx_fatal("g_membed not implemented for openmm");
 127 #endif
 128
 129     t_pargs pa[] = {
 130         { "-xyinit",   FALSE, etREAL,  {&xy_fac},
 131             "Resize factor for the protein in the xy dimension before starting embedding" },
 132         { "-xyend",   FALSE, etREAL,  {&xy_max},
 133             "Final resize factor in the xy dimension" },
 134         { "-zinit",    FALSE, etREAL,  {&z_fac},
 135             "Resize factor for the protein in the z dimension before starting embedding" },
 136         { "-zend",    FALSE, etREAL,  {&z_max},
 137             "Final resize faction in the z dimension" },
 138         { "-nxy",     FALSE,  etINT,  {&it_xy},
 139             "Number of iteration for the xy dimension" },
 140         { "-nz",      FALSE,  etINT,  {&it_z},
 141             "Number of iterations for the z dimension" },
 142         { "-rad",     FALSE, etREAL,  {&probe_rad},
 143             "Probe radius to check for overlap between the group to embed and the membrane"},
 144         { "-pieces",  FALSE,  etINT,  {&pieces},
 145             "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
 146         { "-asymmetry",FALSE, etBOOL,{&bALLOW_ASYMMETRY},
 147             "Allow asymmetric insertion, i.e. the number of lipids removed from the upper and lower leaflet will not be checked." },
 148         { "-ndiff" ,  FALSE, etINT, {&low_up_rm},
 149             "Number of lipids that will additionally be removed from the lower (negative number) or upper (positive number) membrane leaflet." },
 150         { "-maxwarn", FALSE, etINT, {&maxwarn},
 151             "Maximum number of warning allowed" },
 152         { "-start",   FALSE, etBOOL, {&bStart},
 153             "Call mdrun with membed options" },
 154         { "-stepout", FALSE, etINT, {&nstepout},
 155             "HIDDENFrequency of writing the remaining runtime" },
 156         { "-v",       FALSE, etBOOL,{&bVerbose},
 157             "Be loud and noisy" },
 158         { "-mdrun_path", FALSE, etSTR, {&mdrun_path},
 159             "Path to the mdrun executable compiled with this g_membed version" }
 160     };
 161
 162     FILE *data_out;
 163     output_env_t oenv;
 164     char buf[256],buf2[64];
 165     gmx_bool bSucces;
 166
 167     parse_common_args(&argc,argv,0, NFILE,fnm,asize(pa),pa,
 168                       asize(desc),desc,0,NULL, &oenv);
 169
 170     data_out = ffopen(opt2fn("-dat",NFILE,fnm),"w");
 171
 172     fprintf(data_out,"nxy = %d\nnz = %d\nxyinit = %f\nxyend = %f\nzinit = %f\nzend = %f\n"
 173                      "rad = %f\npieces = %d\nasymmetry = %s\nndiff = %d\nmaxwarn = %d\n",
 174                      it_xy,it_z,xy_fac,xy_max,z_fac,z_max,probe_rad,pieces,
 175                      bALLOW_ASYMMETRY ? "yes" : "no",low_up_rm,maxwarn);
 176
 177     fclose(data_out);
 178
 179     sprintf(buf,"%s -s %s -membed %s -o %s -c %s -e %s -nt 1 -cpt -1",
 180                 (mdrun_path==NULL) ? "mdrun" : mdrun_path,
 181                 opt2fn("-f",NFILE,fnm),opt2fn("-dat",NFILE,fnm),opt2fn("-o",NFILE,fnm),
 182                 opt2fn("-c",NFILE,fnm),opt2fn("-e",NFILE,fnm));
 183
 184     if (opt2bSet("-n",NFILE,fnm))
 185     {
 186         sprintf(buf2," -mn %s",opt2fn("-n",NFILE,fnm));
 187         strcat(buf,buf2);
 188     }
 189
 190     if (opt2bSet("-x",NFILE,fnm))
 191     {
 192         sprintf(buf2," -x %s",opt2fn("-x",NFILE,fnm));
 193         strcat(buf,buf2);
 194     }
 195
 196     if (opt2bSet("-p",NFILE,fnm))
 197     {
 198         sprintf(buf2," -mp %s",opt2fn("-p",NFILE,fnm));
 199         strcat(buf,buf2);
 200     }
 201
 202     if (bVerbose)
 203     {
 204         sprintf(buf2," -v -stepout %d",nstepout);
 205         strcat(buf,buf2);
 206     }
 207
 208     if (bStart)
 209     {
 210         printf("Start run with:\n%s\n",buf);
 211         bSucces = system(buf);
 212     }
 213     else
 214     {
 215         printf("You can membed your protein now by:\n%s\n",buf);
 216     }
 217
 218     fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
 219
 220     return 0;
 221 }