src/tools/gmx_membed.c

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