#LyX 1.5.0rc2 created this file. For more info see http://www.lyx.org/ \lyxformat 276 \begin_document \begin_header \textclass revtex4 \begin_preamble \usepackage{listings} \usepackage[]{hyperref} \end_preamble \options rmp, superscriptaddress \language english \inputencoding auto \font_roman times \font_sans helvet \font_typewriter courier \font_default_family default \font_sc false \font_osf false \font_sf_scale 100 \font_tt_scale 100 \graphics default \paperfontsize 10 \spacing single \papersize a4paper \use_geometry false \use_amsmath 1 \use_esint 0 \cite_engine natbib_numerical \use_bibtopic false \paperorientation portrait \secnumdepth 3 \tocdepth 3 \paragraph_separation indent \defskip medskip \quotes_language english \papercolumns 1 \papersides 2 \paperpagestyle default \tracking_changes false \output_changes false \author "Jan H. Meinke" \end_header \begin_body \begin_layout Title \size huge SMMP v. 3 \end_layout \begin_layout Author F. Eisenmenger \end_layout \begin_layout Author Email eisenmenger@fmp-berlin.de \end_layout \begin_layout Affiliation Leibniz-Institut f \begin_inset ERT status inlined \begin_layout Standard { \backslash "u} \end_layout \end_inset r Molekulare Pharmakologie, 13125 Berlin, Germany \end_layout \begin_layout Author U.H.E. Hansmann \end_layout \begin_layout Author Email hansmann@mtu.edu \end_layout \begin_layout Affiliation Michigan Technological University, Houghton, MI 49931-1291, USA \end_layout \begin_layout Affiliation Computational Biology and Biophysics, NIC, 52428 J \begin_inset ERT status inlined \begin_layout Standard { \backslash "u} \end_layout \end_inset lich, Germany \end_layout \begin_layout Author Sh. Hayryan \end_layout \begin_layout Author Email shura@phys.sinica.edu.tw \end_layout \begin_layout Author C.-K. Hu \end_layout \begin_layout Author Email huck@phys.sinica.edu.tw \end_layout \begin_layout Affiliation Institute of Physics, Academia Sinica, Nankang, Taipei 11529, Taiwan \end_layout \begin_layout Author J. H. Meinke \end_layout \begin_layout Author Email j.meinke@fz-juelich.de \end_layout \begin_layout Affiliation Computational Biology and Biophysics, NIC, 52428 J \begin_inset ERT status inlined \begin_layout Standard { \backslash "u} \end_layout \end_inset lich, Germany \end_layout \begin_layout Author S. Mohanty \end_layout \begin_layout Author Email s.mohanty@fz-juelich.de \end_layout \begin_layout Affiliation Computational Biology and Biophysics, NIC, 52428 J \begin_inset ERT status inlined \begin_layout Standard { \backslash "u} \end_layout \end_inset lich, Germany \end_layout \begin_layout Standard \noindent \begin_inset LatexCommand tableofcontents \end_inset \end_layout \begin_layout Section Introduction \end_layout \begin_layout Standard The SMMP ( \series bold S \series default imple \series bold M \series default olecular \series bold M \series default echanics for \series bold P \series default roteins) package is designed for molecular simulations of linear proteins within the standard geometry model. The package contains various modern Monte Carlo algorithms and energy minimizat ion routines. The energy of a protein can be calculated by exploiting one of four force fields: ECEPP/2, ECEPP/3 \begin_inset LatexCommand citep key "ECEPP3" \end_inset , FLEX or Lund. Two subroutines for approximating protein-solvent interactions by means of calculating the solvent accessible surface area of atomic groups are included. Nine sets of solvation parameters are available. All calculations are done with double precision, but the user can easily modify this option by changing the corresponding IMPLICIT statement. \end_layout \begin_layout Standard This manual provides the user with the necessary information for installation of the program, the format of the input and output files, a brief description of the important routines, and describes some limitations of the package. A few example runs providing detailed explanation of input and output are added. More information on SMMP, definitions, background and references for force fields and algorithms as used in the program package can be found in \begin_inset LatexCommand citep key "Eisenmenger2001,Eisenmenger2006" \end_inset . \end_layout \begin_layout Standard SMMP is offered as Open Source. The program is written in standard FORTRAN. We encourage users to improve the code and add their own modules. Please send us your modifications and additions, so we can add them to future releases. Or send us a note and let us know where users can download your modules. The usage and distribution of the program are governed by the following license agreement: \newline \end_layout \begin_layout Quote SMMP is free software; it can be used, modified, and redistributed under the terms of the GNU General Public License version 2 ( \shape italic http://www.gnu.org/licenses/gpl.html \shape default ) as published by the Free Software Foundation, EXCEPT WHERE RESTRICTED BY ONE OF THE FOLLOWING CLAUSES. \newline SMMP shall not be used for commercial or military purposes. The authors request notification if parts of SMMP are incorporated in other programs. \newline The use of SMMP has to be acknowledged in all resulting publications by quoting: \shape italic F. Eisenmenger, U.H.E. Hansmann, S. Hayryan and C.K. Hu, \shape default \series bold [SMMP] A Modern Package for Protein Simulations \series default , \shape italic Comp. Phys. Comm. 138 (2001) 192-212 \shape default ; \series bold An Enhanced Version of SMMP - Open-Source Software Package for Simulation of Proteins \series default , \shape italic Comp. Phys. Comm. 174 (2006) 422- \shape default 429. \emph on J. H. Meinke, S. Mohanty, F. Eisenmenger and U.H.E. Hansmann, \series bold \emph default SMMP v. 3.0 --- Simulating proteins and protein interactions in Python and Fortran \series default , submitted to \emph on Comp.\InsetSpace ~ Phys.\InsetSpace ~ Comm. \emph default \newline SMMP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. \end_layout \begin_layout Standard \noindent The most recent version of SMMP can be found on the program package's web sites: \series bold http://www.smmp05.net \series default or \series bold http://www.phy.mtu.edu/biophys/smmp.htm \series default . Any suggestions for improvements of the code or reports on bugs are welcome. Please send your remarks to the authors. \end_layout \begin_layout Standard \noindent \series bold Acknowledgments: \series default \newline U.H. acknowledges support by a research grant (CHE-9981874) of the National Science Foundation (USA). S.H. and C.K.H. are supported by Academia Sinica (Taipei) and by the National Science Council of the Republic of China (Taiwan) under Contract No. NSC 89-2112-M-001-084. \begin_inset ERT status collapsed \begin_layout Standard \backslash newpage \end_layout \begin_layout Standard \end_layout \end_inset \end_layout \begin_layout Section Changes from Previous Version \end_layout \begin_layout Subsection New Features \end_layout \begin_layout Subsubsection New Force Field \end_layout \begin_layout Paragraph Lund Force Field \end_layout \begin_layout Standard The Lund force fields is a fast force field that uses a minimal set of parameter s to obtain good folding statistics. It has been used successfully for protein folding, unfolding, and aggregation simulations. For details see \begin_inset LatexCommand citep key "Irback2005" \end_inset . \end_layout \begin_layout Subsubsection Multi-molecule simulations \end_layout \begin_layout Standard SMMP 3.0 implements inter-molecular interactions based on the ECEPP/3 force field. \end_layout \begin_layout Subsubsection Biased Gaussian Step \end_layout \begin_layout Standard A new Monte Carlo move that interpolates smoothly between two nearly fixed end points. \end_layout \begin_layout Subsubsection Parallel energy calculation \end_layout \begin_layout Standard The energy calculation for the ECEPP/3 force field can now be done in parallel making it possible to speed up simulations for serial algorithms. Combining the parallel energy calculation with parallel tempering, SMMP can be used efficiently on 2048 processors on an IBM BlueGene/L super computer. \end_layout \begin_layout Subsubsection Python bindings \end_layout \begin_layout Standard We now provide Python bindings to SMMP that make the data structure and algorithms available from (interactive) Python scripts. \end_layout \begin_layout Subsection API Changes \end_layout \begin_layout Standard In previous versions of SMMP the parameters for each algorithm were given in the source code of the algorithm using PARAMETER statements. Now, the parameters are passed as function arguments making the program much more flexible. The following functions and subroutines now have a different signature: anneal, canon, minim, mulcan_par, mulcan_sim, outpdb, outvar, partem_p, and regul. For details see Section\InsetSpace ~ \begin_inset LatexCommand ref reference "sec:Frequently-used-Functions" \end_inset . \end_layout \begin_layout Section Installation \end_layout \begin_layout Standard SMMP requires a FORTRAN 90 compatbile compiler. We have compiled it using gfortran, pgif, xlf, and ifort on various platforms. It should compile and link on any platform with a contemporary FORTRAN compiler. All subroutines are stored in *.f or *.f90 files.An example Makefile is included in the package. There are no machine dependent routines included in SMMP. All COMMON blocks and the limiting parameters are gathered in INCL.H, which is attached to the necessary modules through an INCLUDE statement. After uncompressing and unpacking (via 'tar -xvf' command) the SMMP package into a separate directory For the installation of SMMP the user needs to run \family typewriter make \family default in the directory with the source code. You should edit the Makefile to customize the compiler command and compiler options used. \end_layout \begin_layout Subsection Building pySMMP \end_layout \begin_layout Standard This version of SMMP includes Python bindings to SMMP. To build them, you need f2py, which is part of the numpy package. If f2py is already installed on your system, pySMMP is built by issuing the command \emph on make pybind \emph default . \end_layout \begin_layout Standard \begin_inset ERT status collapsed \begin_layout Standard \backslash newpage \end_layout \begin_layout Standard \end_layout \begin_layout Standard \end_layout \end_inset \end_layout \begin_layout Section \emph on \noun on Biased Gaussian Steps \end_layout \begin_layout Subsection Introduction \end_layout \begin_layout Standard The Bias Gaussian Step, BGS for short, is a semi-local conformational update, involving a concerted rotation of the backbone angles of 4 consecutive residues in the protein chain. The upstream part of the protein backbone, \begin_inset Quotes eld \end_inset before \begin_inset Quotes erd \end_inset the first of the 4 residues used for BGS is left unchanged. The downstream part, \begin_inset Quotes eld \end_inset after \begin_inset Quotes erd \end_inset those 4 residues is moved slightly as a rigid body. The section containing the 4 BGS residues is locally deformed. \end_layout \begin_layout Subsection Parameters \end_layout \begin_layout Standard The increments to the backbone angles involved in the update are first chosen from a Gaussian distribution, and then biased to impart minimal changes to the downstream parts of the chain. If this later \begin_inset Quotes eld \end_inset locality \begin_inset Quotes erd \end_inset condition is relaxed, they would be trivially Gaussian-distributed. There are two parameters of relevance, called \begin_inset ERT status open \begin_layout Standard $a_{ \backslash mathrm{BGS}}$ \end_layout \end_inset and \begin_inset ERT status open \begin_layout Standard $b_{ \backslash mathrm{BGS}}$ \end_layout \end_inset , declared in the common block \noun on bgs_r \noun default . \begin_inset ERT status open \begin_layout Standard $a_{ \backslash mathrm{BGS}}$ \end_layout \end_inset controls the size of the increments, and hence directly influences the acceptance rate of the updates. Smaller the value of \begin_inset ERT status open \begin_layout Standard $a_{ \backslash mathrm{BGS}}$ \end_layout \end_inset , smaller the increments, and hence greater the chances of the update being accepted. Since a trivial zero size update is of little practical interest, even if it is always accepted, it is desirable to have the value of \begin_inset ERT status open \begin_layout Standard $a_{ \backslash mathrm{BGS}}$ \end_layout \end_inset as large as possible. \begin_inset ERT status open \begin_layout Standard $b_{ \backslash mathrm{BGS}}$ \end_layout \end_inset controls the degree of locality of the update, i.e., to what extent the downstrea m parts should be allowed to move. Larger \begin_inset ERT status open \begin_layout Standard $b_{ \backslash mathrm{BGS}}$ \end_layout \end_inset means a stricter locality criterion, and smaller acceptance. These two parameters should not be arbitrarily changed. The default values in SMMP, \begin_inset ERT status open \begin_layout Standard $a_{ \backslash mathrm{BGS}}=300$ \end_layout \end_inset and \begin_inset ERT status open \begin_layout Standard $b_{ \backslash mathrm{BGS}}=10$ \end_layout \end_inset are tuned, to ensure maximum acceptance at a reasonable update size and locality constraint. \end_layout \begin_layout Subsection Detailed balance and BGS \end_layout \begin_layout Standard Information about the current state of the protein is used inside BGS to bias the increments towards greater locality. Thus the current conformation affects the distribution of the incremental angles. This means that the reverse update would have a different probability of being proposed. Because of this, in order to ensure detail balance, an extra multiplicative weight factor is needed for BGS along with the standard Metropolis weight, to decide whether an update should be accepted. The implementation of BGS in SMMP calculates this weight and evaluates whether or not the proposed update was accepted. \end_layout \begin_layout Subsection The function BGS in SMMP \end_layout \begin_layout Standard BGS is implemented in SMMP as a function with the following signature: \end_layout \begin_layout Standard integer function bgs(eol1,dummy) \end_layout \begin_layout Standard The parameter eol1 passed to BGS is the energy of the system before the call. The argument \begin_inset Quotes eld \end_inset dummy \begin_inset Quotes erd \end_inset is a function that calculates the exponent of the Metropolis weight function. The return value is 1 for an accepted BGS move and 0 for a rejected move. \end_layout \begin_layout Subsubsection Using BGS in SMMP \end_layout \begin_layout Standard To use BGS, one must ensure that the subroutine init_lund, which initialises several arrays used in BGS, is called during the initialisation phase of the program. These arrays store information about the BGS-capable angles in the protein chain. Therefore, init_lund should be called again, in case the program changes the amino acid sequence of a molecule or the total number of molecules in the system. \end_layout \begin_layout Standard The global switch \noun on upchswitch \noun default controls how BGS is used by the metropolis subroutine, which is normally the routine that calls BGS. For \noun on upchswitch=0, BGS \noun default is never used. For \noun on upchswitch=1, \noun default it is used with a probabilty \noun on bgsprob \noun default , which is another global parameter. For \noun on upchswitch=2, \noun default BGS is used when possible with a temperature dependent probability. At higher temperatures abrupt single angle changes are more effective in browsing the conformation space quickly. A local update like BGS should therefore have a small probability at high temperatures. At low temperatures, often the system is in a compact state with significant secondary structure, and it is desirable to explore the neighbourhood of that state to refine the structure. BGS is more efficient at this compared to single angle updates, and hence should have a higher probability at lower temperatures. The choice, \noun on upchswitch=2 \noun default uses this scheme for BGS use. \end_layout \begin_layout Section Getting Started \end_layout \begin_layout Standard The quickest way to get started is to use one of the examples in the EXAMPLES directory as a template. The examples are described in Section\InsetSpace ~ \begin_inset LatexCommand ref reference "sec:Examples" \end_inset . \end_layout \begin_layout Standard SMMP does not include an interpretor of user-defined commands instead the preparation of a simulation is done in the MAIN module calling corresponding simulation subroutine(s). After adjusting MAIN to your needs, SMMP has to be re-compiled. We provide a sample main file with detailed comments. For a comment-free version see Listing \begin_inset LatexCommand ref reference "listing:mainmc" \end_inset . \end_layout \begin_layout Standard \begin_inset listings lstparams "float=h,language=Fortran,numbers=left" inline false status open \begin_layout Standard program main \end_layout \begin_layout Standard include "../INCL.H" \end_layout \begin_layout Standard \end_layout \begin_layout Standard character*80 libdir, seqfile, varfile \end_layout \begin_layout Standard character grpn*4,grpc*4 \end_layout \begin_layout Standard integer nequi, nsweeps, nmes \end_layout \begin_layout Standard double precision tmax, tmin \end_layout \begin_layout Standard logical lrand \end_layout \begin_layout Standard \end_layout \begin_layout Standard libdir = '../SMMP/' \end_layout \begin_layout Standard ientyp = 0 \end_layout \begin_layout Standard sh2 = .true. \end_layout \begin_layout Standard itysol = 0 \end_layout \begin_layout Standard call init_energy(libdir) \end_layout \begin_layout Standard \end_layout \begin_layout Standard grpn = 'nh2' \end_layout \begin_layout Standard grpc = 'cooh' \end_layout \begin_layout Standard iabin = 1 \end_layout \begin_layout Standard seqfile='EXAMPLES/enkefa.seq' \end_layout \begin_layout Standard varfile='EXAMPLES/enkefa.var' \end_layout \begin_layout Standard ntlml = 0 \end_layout \begin_layout Standard call init_molecule(iabin,grpn,grpc,seqfile,varfile) \end_layout \begin_layout Standard \end_layout \begin_layout Standard seed = 81236 \end_layout \begin_layout Standard call sgrnd(seed) \end_layout \begin_layout Standard upchswitch= 0 \end_layout \begin_layout Standard \end_layout \begin_layout Standard nequi=100 \end_layout \begin_layout Standard nsweeps=100000 \end_layout \begin_layout Standard nmes=1000 \end_layout \begin_layout Standard temp = 300 \end_layout \begin_layout Standard lrand=.true. \end_layout \begin_layout Standard call canon(nequi, nsweeps, nmes, temp,lrand) \end_layout \begin_layout Standard \end_layout \begin_layout Standard end program main \end_layout \begin_layout Standard \begin_inset Caption \begin_layout Standard The basic main function for a Monte Carlo Simulation \begin_inset LatexCommand label name "listing:mainmc" \end_inset \end_layout \end_inset \end_layout \end_inset \end_layout \begin_layout Standard The following steps summarize how to run SMMP: \end_layout \begin_layout Enumerate Assign the path to the directory containing the standard amino acid residue libraries and the file 'charges' to the character variable 'libdir' (Line 10). \end_layout \begin_layout Enumerate Select the force field and solvation model by setting the four paramters 'flex', 'sh2', 'epsd', and 'itysol' to their appropriate values: \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard * ientyp = 0 : ECEPP \end_layout \begin_layout Standard 1 : FLEX \end_layout \begin_layout Standard 2 : LUND \end_layout \begin_layout Standard 3 : ECEPP + Abagyan \end_layout \begin_layout Standard * sh2 =.TRUE. : ECEPP/2 potential, sh2=.FALSE.: ECEPP/3 \end_layout \begin_layout Standard * epsd=.TRUE. : Distant dependent dielectric permittivity \end_layout \begin_layout Standard epsd=.FALSE. : epsilon = 2 \end_layout \begin_layout Standard * itysol = 0 : gas phase, \end_layout \begin_layout Standard itysol > 0 : approximation of protein-solvent interactions by means \end_layout \begin_layout Standard of a solvent accessible surface area approach with \end_layout \begin_layout Standard numerical estimation of the accessible area, \end_layout \begin_layout Standard itysol < 0 : same as above, but the accessible area is calculated \end_layout \begin_layout Standard analytically (considerably than itysol > 0). \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset (Line 11--13). \end_layout \begin_layout Enumerate \shape italic call init_energy(libdir) \shape default to initialize the energy parametrization (Line 12). \end_layout \begin_layout Enumerate Choose the N-terminal and C-terminal groups by assigning to appropriate strings to variables 'grpn' and 'grpc' (Line 16--17) \end_layout \begin_layout Enumerate Choose how the initial input is read in: \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard * iabin = 0 : read from PDB-file \end_layout \begin_layout Standard * iabin = 1 : read from sequence (and configuration) file \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset (Line 18) \end_layout \begin_layout Enumerate Enter the names of the corresponding file(s) ' \shape italic seqfil \shape default ' and ' \shape italic varfil \shape default ' (Line 19--20). If there are empty, they'll be requested through an interactive command line dialog, issued by ' \shape italic call init_molecule \shape default ' (Line 22). The user can easily just supply the corresponding file names in a "here document", as shown in the example scripts . \end_layout \begin_layout Enumerate At this point the program is ready for calling task subroutines. Here we perform a canonical Monte Carlo simulation at T=300K. Two kinds of updates are avaible: a single angle update and the biased Gaussian step (bgs). For this simulation we turn bgs off by setting \family typewriter upchswitch=0 \family default . Detailed examples for this and other simulation tasks can be found in the following section. Normally, simulation subroutines write data in output files, but one can also use output routines such as ` \shape italic outpdb \shape default ' in ' \shape italic main \shape default '. The minimal output (written into standard output) is the name of the sequence file (extension .seq), name of configuration file (extension .var), and for each residue a list of dihedral angles together with their initially assigned values. \end_layout \begin_layout Enumerate For parallel tempering jobs on on a multiprocessor system one has to replace ' \shape italic main \shape default ' by ' \shape italic pmain \shape default '. The above protocol still applies, but one has to provide in addition the number of nodes with the variable ' \shape italic no \shape default '. \end_layout \begin_layout Section Programming with SMMP \end_layout \begin_layout Standard While SMMP provides several efficient sampling algorithms, you may want to improve an existing or implement a new algorithm. This section goes into more details and explains some of the nuts and bolts of SMMP. \end_layout \begin_layout Paragraph Where can I find information about atoms, residues, proteins, and dihedrals \end_layout \begin_layout Standard Configurations are stored in one dimensional arrays defined in INCL.H and shared through common blocks. \begin_inset Note Note status collapsed \begin_layout Standard I hope, we'll have moved most of the arrays to modules before publication. \end_layout \end_inset It helps to think somewhat object oriented to understand the structure of the arrays. In general integer arrays start with the letter 'i'. The abbreviation 'ml' relates to molecules, 'rs' to residues, and 'at' to atoms. Several boundary arrays are pointers to the first and last element and end in '1' for the first element and '2' for the last element. For example, irsml1(i) is the index of the first residue in molecule \begin_inset Formula $i$ \end_inset , iatrs1(irsml(i)) is the index of the first atom of the first residue of molecule \begin_inset Formula $i$ \end_inset . If iat = iatrs1(irsml(i)), then we can access the coordinates of the atom using xat(iat), yat(iat), z(iat). Using this scheme we can cycle through all the atoms in the system. An example is shown in Listing \begin_inset ERT status open \begin_layout Standard \backslash ref{listing:outpdb} \end_layout \end_inset . \begin_inset ERT status open \begin_layout Standard \backslash lstset{language=Fortran} \end_layout \begin_layout Standard \backslash begin{lstlisting}[float=p, caption={Using the indices to cycle through all the atoms in the system. This snippet is taken from the subroutine outpdb, which writes the current configuration to a file using the pdb format.}, label=listing:outpdb] \end_layout \begin_layout Standard irs=0 \end_layout \begin_layout Standard iat0 = 0 \end_layout \begin_layout Standard do inml = nml0, nml1 \end_layout \begin_layout Standard iat = iat0 \end_layout \begin_layout Standard chid = char(64 + inml) \end_layout \begin_layout Standard ifirs=irsml1(inml) \end_layout \begin_layout Standard ifiat=iatrs1(ifirs) \end_layout \begin_layout Standard do nrs=ifirs,irsml2(inml) \end_layout \begin_layout Standard irs=irs+1 \end_layout \begin_layout Standard res(1:)=seq(nrs)(1:3) \end_layout \begin_layout Standard do i=iatrs1(nrs),iatrs2(nrs) \end_layout \begin_layout Standard iat=iat+1 \end_layout \begin_layout Standard atnm=' ' \end_layout \begin_layout Standard atnm(2:5)=nmat(i) ! nmat(i) contains the name of atom i \end_layout \begin_layout Standard if (atnm(5:5).ne.' ') then \end_layout \begin_layout Standard atnm(1:1)=atnm(5:5) !! \end_layout \begin_layout Standard atnm(5:5)=' ' \end_layout \begin_layout Standard endif \end_layout \begin_layout Standard z=zat(i)-rz \end_layout \begin_layout Standard call toupst(atnm) ! converts a string to all uppercase \end_layout \begin_layout Standard call toupst(res) \end_layout \begin_layout Standard \end_layout \begin_layout Standard write (npdb,1) 'ATOM',iat,atnm,res(1:3),chid,irs,cdin, \end_layout \begin_layout Standard xat(i), yat(i),zat(i),occ,bva \end_layout \begin_layout Standard enddo \end_layout \begin_layout Standard enddo \end_layout \begin_layout Standard enddo \end_layout \begin_layout Standard \backslash end{lstlisting} \end_layout \end_inset \begin_inset Float table placement p wide false sideways false status collapsed \begin_layout Standard \begin_inset Tabular \begin_inset Text \begin_layout Standard Variable Name \end_layout \end_inset \begin_inset Text \begin_layout Standard Type \end_layout \end_inset \begin_inset Text \begin_layout Standard Description \end_layout \end_inset \begin_inset Text \begin_layout Standard ityat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of int \end_layout \end_inset \begin_inset Text \begin_layout Standard Index of type of atom \end_layout \end_inset \begin_inset Text \begin_layout Standard xat, yat, zat \end_layout \end_inset \begin_inset Text \begin_layout Standard arrays of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Cartesian coordinates of atom \end_layout \end_inset \begin_inset Text \begin_layout Standard cgat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Partial electric charge on atom \end_layout \end_inset \begin_inset Text \begin_layout Standard nmat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of character*4 \end_layout \end_inset \begin_inset Text \begin_layout Standard Name of atom, e.g., C, N, O, H1, H2 \end_layout \end_inset \begin_inset Text \begin_layout Standard iowat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of int \end_layout \end_inset \begin_inset Text \begin_layout Standard Index of preceeding atom (see Fig.\InsetSpace ~ \begin_inset LatexCommand ref reference "fig:Enumeration-of-atoms" \end_inset ), 0 if no preceeding atom. \end_layout \end_inset \begin_inset Text \begin_layout Standard iyowat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of int \end_layout \end_inset \begin_inset Text \begin_layout Standard Type of bond with preceeding atom. \begin_inset Note Note status open \begin_layout Standard Maybe, this should be called ityinbdat. \end_layout \end_inset \end_layout \end_inset \begin_inset Text \begin_layout Standard nbdat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of int \end_layout \end_inset \begin_inset Text \begin_layout Standard Number of bonds to succeeding atom (outgoing bonds) \end_layout \end_inset \begin_inset Text \begin_layout Standard ibdat \end_layout \end_inset \begin_inset Text \begin_layout Standard two dimensional array of int \end_layout \end_inset \begin_inset Text \begin_layout Standard Indices of successor atoms \end_layout \end_inset \end_inset \end_layout \begin_layout Standard \begin_inset Caption \begin_layout Standard \begin_inset LatexCommand label name "tab:atomic-properties" \end_inset List of variables storing atomic properties. All the properties of atom \begin_inset Formula $i$ \end_inset are accessed as the \begin_inset Formula $i$ \end_inset th element of the arrays. \end_layout \end_inset \end_layout \end_inset \begin_inset Float table placement p wide false sideways false status collapsed \begin_layout Standard \begin_inset Tabular \begin_inset Text \begin_layout Standard Name \end_layout \end_inset \begin_inset Text \begin_layout Standard Type \end_layout \end_inset \begin_inset Text \begin_layout Standard Description \end_layout \end_inset \begin_inset Text \begin_layout Standard xbaat, ybaat, zbaat \end_layout \end_inset \begin_inset Text \begin_layout Standard 3 arrays of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Rotation axis for valence angle. (see Fig. ) \end_layout \end_inset \begin_inset Text \begin_layout Standard baat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Valence angle in radians \end_layout \end_inset \begin_inset Text \begin_layout Standard xtoat, ytoat, ztoat \end_layout \end_inset \begin_inset Text \begin_layout Standard 3 arrays of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Rotation axis for dihedral (torsion) angle. (see Fig.) \end_layout \end_inset \begin_inset Text \begin_layout Standard toat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Torsion angle defined by iowat(iowat(iowat(i))), iowat(iowat(i)), iowat(i), i in radians (see Fig.) \end_layout \end_inset \begin_inset Text \begin_layout Standard blat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of double \end_layout \end_inset \begin_inset Text \begin_layout Standard Length of bond. \end_layout \end_inset \begin_inset Text \begin_layout Standard ixmsat \end_layout \end_inset \begin_inset Text \begin_layout Standard array of int \end_layout \end_inset \begin_inset Text \begin_layout Standard Index of moving set belonging to atom \end_layout \end_inset \begin_inset Text \begin_layout Standard \end_layout \end_inset \begin_inset Text \begin_layout Standard \end_layout \end_inset \begin_inset Text \begin_layout Standard \end_layout \end_inset \end_inset \end_layout \begin_layout Standard \begin_inset Caption \begin_layout Standard \begin_inset LatexCommand label name "tab:dihedral-to-atom" \end_inset Each dihedral and valence angle and the corresponding moving set is associated with a covalent bond which in turn is defined by its tail atom. This table list the related variables. \end_layout \end_inset \end_layout \end_inset \begin_inset Float figure wide false sideways false status collapsed \begin_layout Standard \align center \begin_inset Graphics filename angle_defs.eps scaleBeforeRotation subcaption \end_inset \begin_inset Graphics filename dihedral_defs.eps scaleBeforeRotation subcaption \end_inset \end_layout \begin_layout Standard \begin_inset Caption \begin_layout Standard \begin_inset LatexCommand label name "fig:angle-defs-in-SMMP" \end_inset Angle in SMMP. There are two different types of angles in SMMP (and proteins in general): bond and torsional angles. The bond or valence angle is show in (a). The red arrow indicates the rotation axis and is proportional to the cross product of the two bonds \begin_inset Formula $(1,4)$ \end_inset and \begin_inset Formula $(4,25)$ \end_inset making up the angle. The angle \begin_inset Formula $\alpha$ \end_inset is associated with atom \begin_inset Formula $4$ \end_inset . The dihedral or torsional angles are rotations around a bond and determined by the angles between the two adjacent planes. For example, the angle \begin_inset Formula $\psi$ \end_inset in (b) is the angle between the plane defined by atoms 1, 4, and 25, and the plane defined by 4, 25, and 28. We could therefore reference this angle using the quadruplet (1, 4, 25, 28). Note that atom 1 is the predecessor of atom 4, atom 4 the predecessor of 25, and, finally, atom 25 the predecessor of 28. SMMP uses this relation to associate the angle \begin_inset Formula $\psi$ \end_inset with atom 28. \end_layout \end_inset \end_layout \end_inset \end_layout \begin_layout Standard Table \begin_inset LatexCommand ref reference "tab:atomic-properties" \end_inset lists the arrays containing the atomic properties such as charge, position, and in- and outgoing bonds. In the standard geometry configurations can defined by the dihedral angles. The variables associated with the angles are given in Table \begin_inset LatexCommand ref reference "tab:dihedral-to-atom" \end_inset and a description of the angles can be found in Fig.\InsetSpace ~ \begin_inset LatexCommand ref reference "fig:angle-defs-in-SMMP" \end_inset . \end_layout \begin_layout Paragraph \begin_inset Note Note status open \begin_layout Standard Why are the names so short and what do they mean? \end_layout \begin_layout Paragraph Setting up a force field \end_layout \begin_layout Paragraph Loading a molecule from a .seq and .var file \end_layout \begin_layout Standard fixed variables, default values, access? \end_layout \begin_layout Paragraph Loading a molecule from a PDB file \end_layout \begin_layout Standard where are things stored? what does regularization do? \end_layout \end_inset \end_layout \begin_layout Section Limitations \end_layout \begin_layout Standard All parameters which limit the usage of SMMP are stored in the file INCL.H. The most important ones are listed below. The values of these parameters should be changed only in a consistent way! \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxml=1 max. number of molecules \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxrs=100 max. total number of residues \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxat=2000 max. total number of atoms \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxbd=3 max. number of bonds to following atoms \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxvr=mxrs*5 max. number of local variables \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxms=mxvr*3 max. total number of moving sets \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxvw=mxat*4 max. number of vdw domains \end_layout \begin_layout Standard \end_layout \begin_layout Standard mx14=mxat*4 max. number of '1-4' partners \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxath=100, max. number of atoms in help-arrays \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxvrh=mxath max. number of variables in help \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxtyat=18 max. number of energetic atom-types \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxhbdo=4 max. types of Hydrogens as donors in HB \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxhbac=6 max. types of atoms as acceptors in HB \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxtyto=19 max. number of types of torsional potentials \end_layout \begin_layout Standard \end_layout \begin_layout Standard nrsty=35 max. number of residue types \end_layout \begin_layout Standard \end_layout \begin_layout Standard mxtysol=9 the number of solvation parameters sets \end_layout \begin_layout Standard \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard \noindent Note also the following restrictions in the current version of SMMP: \end_layout \begin_layout Enumerate A single amino acid residue can not be simulated with the FLEX potential \end_layout \begin_layout Enumerate A protein must not start with a prolyl residue. \end_layout \begin_layout Section Libraries and Parameter Files \end_layout \begin_layout Standard The residues which can be used with each parameter set are described in files \shape italic lib.sh2, lib.sh3 \shape default , and \shape italic lib.flex \shape default , respectively. The file \shape italic charges \shape default is needed for N- and C-terminal residues with FLEX parameters. The name of the directory, containing these 4 files should be given in string 'libdir', which is assigned in module MAIN. The choice between parameter sets for potentials ECEPP/2 and ECEPP/3 is done by the logical variable ' \shape italic sh2 \shape default '. The potential FLEX is set by the logical variable ' \shape italic flex \shape default '. \end_layout \begin_layout Standard The amino acid residue libraries contain chemical and structural data for each residue according the IUPAC-IUB nomenclature. The library files consist of blocks of records with each block representing one residue. The first line in each block starts with a '#' and contains the name of the residue and its total number of atoms. Each following line within the block describes one atom. Listed are: its name, the bond length and the valence angle to construct this atom, the type and the name of corresponding torsion angle, the partial charge for the atom, the type of the atom, as well as the previous and the following atoms (max. 3), it is connected to. The atom and torsional parameters in the libraries are the same as in the original ECEPP/2/3 and FLEX potentials. However, it should be noted that SMMP has its own numeration of the atom types and the types of torsions (see [1,2] for more details). \end_layout \begin_layout Standard The arrays with the parameters for non bonded pairwise potentials (Van der Waals, hydrogen bonds, electrostatic) and the atomic solvation parameters are stored in the BLOCK DATA section of ' \shape italic init_energy.f \shape default '. \end_layout \begin_layout Standard SMMP employs 9 different sets of atomic solvation parameters. We do not include references to the original works where these sets are defined in this document, since this manual is only a supplement to our articles \begin_inset LatexCommand cite key "Eisenmenger2001,Eisenmenger2006" \end_inset where all relevant references can be found. The set of solvation parameters chosen is indicated by the value of integer variable ' \shape italic itysol \shape default ' in the MAIN module. If setting \shape italic itysol \shape default to a positive value between \shape italic 1 \shape default and \shape italic 9 \shape default in MAIN, the numerical method \shape italic enysol \shape default is applied to compute the solvation energy term. By choosing a negative integer value for \shape italic itysol \shape default from the interval \shape italic [-1,-9] \shape default , the analytical method \shape italic esolan \shape default is used to calculate the solvation energy. The tesselation points for numerically calculating the solvent accessible surface area have to be provided in an external file (named \shape italic tes.dat \shape default in our package) and read in during initialization of solvation parameters in \shape italic init_energy.f \shape default . \end_layout \begin_layout Section Input Files \end_layout \begin_layout Standard SMMP requires an input file that specifies the sequence of amino acid residues of the protein in plain ASCII-format. \end_layout \begin_layout Subsection PDB-File \begin_inset LatexCommand label name "sub:PDB-File" \end_inset \end_layout \begin_layout Standard Setting the variable ' \shape italic iabin = 0 \shape default ' in \shape italic main.f \shape default the sequence of amino acids is read from a PDB-file - the standard format in which protein structures are deposited in the Protein Data Bank. The atomic coordinates are also read from the PDB-file and can serve as a start configuration. However, the fixed bond lengths and angles in the standard geometry will slightly differ from the actual bond lengths and angles in the PDB-structure. Forcing the molecule into standard bonding geometry corresponding to a given potential may lead to un-physically high energies. Regularization through \shape italic regul.f \shape default allows to obtain an optimized structure with standard bonding geometry with low internal energy without differing significantly from the PDB structure. \end_layout \begin_layout Subsection Sequence File \begin_inset LatexCommand label name "sub:Sequence-File" \end_inset \end_layout \begin_layout Standard Setting ` \shape italic iabin = 1 \shape default ' the amino acid sequence is read from a separate sequence file. Its first line must start with a '#' and may contain the name for the molecule. The residues in the following lines should be named as in the library files \shape italic lib.sh2, lib.sh3 \shape default , or \shape italic lib.flex \shape default . Residue names are not case-sensitive and should be separated from each other by at least one white space. \end_layout \begin_layout Standard Currently, the following residue types can be used in SMMP: \end_layout \begin_layout Enumerate Neutral: \newline ala, arg, asn, asp, cys, gln, glu, gly, his* ; hise* ; hyp*, ile, leu, lys, met, phe, pro* , ser, thr, trp, tyr, val \newline (where the residues marked by an asterix are characterized by: \newline his - hydrogen at N-delta atom \newline hise - hydrogen at N-epsilon atom \newline pro, hyp - down-puckering) \end_layout \begin_layout Enumerate Charged: \newline arg+, asp-, glu-, his+, lys+. \end_layout \begin_layout Enumerate Variants, only with ECEPP/3 parameters: \newline hypu, prou - up-puckering \newline cpro, cpru - cis-Pro with down-, up-puckering \newline pron, pro+ - N-terminal neutral, charged with respect to NH2+ Pro \end_layout \begin_layout Standard End-groups need to be specified in \shape italic main.f \shape default and are added through subroutine \shape italic addend \shape default to the N- and C-terminus, correspondingly. Possible values for the end groups are nh2 and nh3+ for the N terminal and cooh and coo- for the C terminal. Note, that the acetyl group (' \shape italic ace \shape default ') and N'methylamide (' \shape italic nme \shape default ') need to be added as \shape italic residues \shape default in the sequence file. In that case, input into ' \shape italic addend \shape default ' will be ignored by SMMP. Note also, that at present, molecules with only one single residue cannot be simulated with the FLEX data set and the first residue should not be a prolyl residue. \end_layout \begin_layout Standard An example sequence file (provided with SMMP) may look as follows: \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard ------------------- enkefa.seq -------------------------- \end_layout \begin_layout Standard # Met-Enkephalin \end_layout \begin_layout Standard Tyr Gly Gly Phe Met \end_layout \begin_layout Standard --------------------------------------------------------- \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Subsection Configuration File (Dihedrals) \begin_inset LatexCommand label name "sub:Configuration-File-(Dihedrals)" \end_inset \end_layout \begin_layout Standard With option ' \shape italic iabin = 1 \shape default ' SMMP allows to provide a designated start configuration by listing dihedral angles in a second input file. If no name of any configuration file is given (or the name of a non-existing file is entered), all variables retain their default values as defined in the libraries. \end_layout \begin_layout Standard This configuration file has to follow a fixed structure in which each line can be considered as a COMMAND for subroutine ' \shape italic redvar \shape default '. The following syntax has to be used in each COMMAND: RESIDUE : VARIABLE : VALUE Hence, each line can consist of up to 3 FIELDS, separated by an ':'. In the first field the RESIDUE is selected through an INTEGER number which marks its position in the amino acid sequence. The second field should contain a string with the name of the VARIABLE, i.e. names the specific dihedral angle. The last field provides the value for the VARIABLE (a REAL number) and is mandatory. A symbol '&' following this number indicates that the variable will be fixed to that value throughout the whole simulation process. Spaces are not significant and are therefore ignored. Empty COMMANDS, as ' : : ', and empty lines or lines containing '#' are ignored. Missing fields are interpreted as 'for all': \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard 1 : phi : 180 Set phi of residue #1 to 180 degrees \end_layout \begin_layout Standard 1 : phi : 180& Keep phi of residue #1 fixed to 180 degrees \end_layout \begin_layout Standard phi : 180 Set phi of ALL residues to 180 degrees \end_layout \begin_layout Standard 180 Set all variables to 180 degrees \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset Several consecutive residues or variables can be indicated by ZONES of indices. Several NAMES can be indicated by wild-card ('*') and are case-sensitive: \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard 1-4 : phi : 180 Set 'phi' of residues #1,#2,#3, and #4 to 180 deg. \end_layout \begin_layout Standard 5 : x* : 60 Set all xi-angles of residue #5 to 60 deg. \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset Several NAMES or INDICES may be given in the same field, when separated by ','. Similarly, several commands may be given on the same line and must be separated by ';' \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard phi, psi : -30 Set all phi & psi to -30 deg. \end_layout \begin_layout Standard phi:-65; psi:-45 Set all phi=-65, all psi=-45 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset Note, that dihedral angles are defined slightly differently in SMMP than in standard ECEPP (residue index in parenthesis): \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard - psi of the i-th residue = defined by N (i-1) - CA(i-1) - C (i-1) - N(i) \end_layout \begin_layout Standard - omg of the i-th residue = defined by CA(i-1) - C (i-1) - N (i) - CA(i) \end_layout \begin_layout Standard - phi for the i-th residue = defined by C (i-1) - N (i) - CA(i) - C(i) \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset except: \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard phi for the 1st residue = defined by H1(1) - N (1) - CA(1) - C(1) \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset (H2-atom in the N-terminal NH2-group is added via dihedral 'H2-N-CA-C'= phi + 120) \end_layout \begin_layout Standard If a molecule consists of n residues then: \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard - pst = defines dihedral N (n) - CA(n) - C (n) - Oxt(n) \end_layout \begin_layout Standard - omt = defines dihedral CA(n) - C (n) - Oxt(n) - Hxt(n) \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset (Oxt & Hxt comprise the hydroxyl of the C-terminal carboxyl group) \end_layout \begin_layout Standard The following example configuration file \shape italic enkefa.var \shape default is provided with the program and may be used with the sequence file \shape italic enkefa.seq \shape default : \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard ------------------- enkefa.var -------------------------- \end_layout \begin_layout Standard 1 : phi : -86.24 \end_layout \begin_layout Standard 1 : x1 : -172.59 \end_layout \begin_layout Standard 1 : x2 : 78.71 \end_layout \begin_layout Standard 1 : x6 : -165.88 \end_layout \begin_layout Standard 2 : psi : 156.18 \end_layout \begin_layout Standard 2 : omg : 180.00 \end_layout \begin_layout Standard 2 : phi : -154.53 \end_layout \begin_layout Standard 3 : psi : 83.64 \end_layout \begin_layout Standard 3 : omg : 180.00 \end_layout \begin_layout Standard 3 : phi : 83.66 \end_layout \begin_layout Standard 4 : psi : -73.86 \end_layout \begin_layout Standard 4 : omg : 180.00 \end_layout \begin_layout Standard 4 : phi : -137.04 \end_layout \begin_layout Standard 4 : x1 : 58.79 \end_layout \begin_layout Standard 4 : x2 : 94.60 \end_layout \begin_layout Standard 5 : psi : 19.33 \end_layout \begin_layout Standard 5 : omg : -180.00 \end_layout \begin_layout Standard 5 : phi : -163.63 \end_layout \begin_layout Standard 5 : x1 : 52.76 \end_layout \begin_layout Standard 5 : x2 : 175.28 \end_layout \begin_layout Standard 5 : x3 : -179.83 \end_layout \begin_layout Standard 5 : x4 : -58.57 \end_layout \begin_layout Standard 5 : pst : 160.45 \end_layout \begin_layout Standard 5 : omt : 180.00 \end_layout \begin_layout Standard --------------------------------------------------------- \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard \begin_inset Float figure wide false sideways false status open \begin_layout Standard \align center \begin_inset Graphics filename atom_numbering.eps scaleBeforeRotation \end_inset \end_layout \begin_layout Standard \begin_inset Caption \begin_layout Standard \begin_inset LatexCommand label name "fig:Enumeration-of-atoms" \end_inset Enumeration of atoms in SMMP. Atoms are numbered starting from the N terminal of the protein. The terminal nitrogen has the index 1. Counting continues with the two attached hydrogens and then moves along the backbone. After the C \begin_inset Formula $_{\alpha}$ \end_inset , we count along the side chain following a breadth-first scheme. Rings are enumerated counter clockwise. \end_layout \end_inset \end_layout \end_inset \end_layout \begin_layout Section Frequently used Functions and Subroutines \begin_inset LatexCommand label name "sec:Frequently-used-Functions" \end_inset \end_layout \begin_layout Subsection Task Subroutines \end_layout \begin_layout Standard In the following we describe some important simulation subroutines and their requisites. For more information see [1,2]. Following the philosophy behind SMMP these routines are rather simple and may serve as templates which the user can modify according to his special needs. \end_layout \begin_layout Subsubsection Subroutine init_molecule(iabin, grpn, grpc, seqfile, varfile) \end_layout \begin_layout Standard (Usage: ' \emph on call init_molecule(iabin, grpn, grpc, seqfile, varfile)') \end_layout \begin_layout Standard Initializes a new molecule and adds it to the system. The information is either read from a sequence (see \begin_inset LatexCommand ref reference "sub:Sequence-File" \end_inset ) and, optionally, a variable file (see \begin_inset LatexCommand ref reference "sub:Configuration-File-(Dihedrals)" \end_inset ) if \emph on iabin=1 \emph default or from a pdb file (see \begin_inset LatexCommand ref reference "sub:PDB-File" \end_inset ) if \emph on \begin_inset Formula $iabin\neq1$ \end_inset \emph default . If the molecule is read from a pdb file, the name of the file is given in \emph on seqfile \emph default . By repeatedly calling init_molecule, you can add several interacting molecules. The N-terminal group is determined by \emph on grpng \emph default and the C-terminal group \emph on grpc. \end_layout \begin_layout Subsubsection Subroutine REGUL(nml, iter, nsteps, \begin_inset Formula $\varepsilon$ \end_inset ) \begin_inset LatexCommand index name "REGUL" \end_inset (file regul.f) \end_layout \begin_layout Standard (Usage: ' \shape italic call regul(nml, iter, nsteps, \begin_inset Formula $\varepsilon$ \end_inset ) \shape default ') \end_layout \begin_layout Standard For a given PDB-structure the actual bond lengths and bond angles generally differ from those in the standard geometry model assumed with the ECEPP or FLEX potential. This subroutine is applied to determine the configuration within the standard geometry that is as close as possible to a given PDB-structure with an energy as low as possible. This process of fitting a PDB structure into standard geometry is called regularization. ' \shape italic regul \shape default ' requires the PDB-structure to be read into SMMP by setting ' \shape italic iabin = 0 \shape default ' in ' \shape italic main \shape default '. With this option, ` \shape italic init_molecule \shape default ' calls the subroutines ` \shape italic pdbread \shape default ' and ` \shape italic pdbvars \shape default ' that also to measure the dihedral angles in the PDB-structure. Using these dihedral angles a first model of the molecule is built using standard bonding geometry. Hence, ' \shape italic regul \shape default ' has to be called AFTER ' \shape italic init_molecule \shape default '. ' \shape italic nml \shape default ' has to be set to 1 in the present version of SMMP. \end_layout \begin_layout Standard ' \shape italic regul \shape default ' starts with a minimization of the initial structure using only the sum of squared distances between atom positions of the SMMP structure to the PDB-structure as \begin_inset Quotes erd \end_inset constraint energy \begin_inset Quotes erd \end_inset . After this initial step, a list of bad atom contacts, i.e. atoms with vdW-energy of more than 2 kcal/mol (as calculated by ' \shape italic cnteny \shape default ') and the root-mean-square-deviation (rmsd) to the PDB-structure is printed out. In a second step, the physical energy is minimized, allowing only hydrogens to move, and bad contacts and the rmsd are evaluated, again. A number of iterations follow where a composite energy as the weighted sum of physical and constraint energy is minimized. In this iterations the weight of the constraint energy is successively lowered to zero, and that of the physical energy raised to one. At the end, the dihedral angles of the final structure are printed out together with its rmsd and list of remaining bad contacts. \end_layout \begin_layout Subsubsection Subroutine MINIM(imin, maxit, \begin_inset Formula $\varepsilon$ \end_inset ) \begin_inset LatexCommand index name "MINIM" \end_inset (file minim.f): \end_layout \begin_layout Standard (Usage: \shape italic 'call minim(imin \shape default , maxit, \begin_inset Formula $\varepsilon$ \end_inset )') \newline This subroutine minimizes the current configuration of the molecule. It can be called from any point during the simulation once a configuration is defined. When calling this subroutine one needs to specify the type of minimizer through variable ' \shape italic imin \shape default '. Setting ' \shape italic imin = 1 \shape default ' chooses a Quasi-Newton minimizer by calling ' \shape italic minqsn \shape default ', ' \shape italic imin = 2 \shape default ' a Conjugate Gradient minimizer (' \shape italic mincjg \shape default '). The performance of these minimizers can be controlled by setting various parameters in MINIM. The two most commonly used ones, the accuracy (' \begin_inset Formula $\varepsilon$ \end_inset ') and the total number of function calls 'mxit' are passed as arguments. \end_layout \begin_layout Subsubsection Subroutine CANON(nequi, nswp, nmes, temp, lrand) \begin_inset LatexCommand index name "CANON" \end_inset (file canon.f): \end_layout \begin_layout Standard (Usage: ' \emph on call canon(nequi, nswp, nmes, temp, lrand) \emph default ') \newline The subroutine carries out a canonical simulation at temperature 'temp'. The simulation starts from a random configuration if the logical parameter 'lrand' is set to \shape italic .TRUE. \shape default . Otherwise the simulation starts from the current conformation of dihedral angles as stored in array 'vlvr'. One also has to set the number of Monte Carlo sweeps for the simulation 'nswp', and the equilibration 'nequi', as well as the number 'nmes' of Monte Carlo sweeps between measurements. Three output files are created: files `start.pdb' and 'final.pdb' are for storing the start and final configurations of the molecule, respectively, in a data format compatible with Protein Data Bank. The file 'time.d' contains the time series of some quantities. In our template we measure and store in this file the total energy 'eol'; the radius-of-gyration `rgy' (by means of subroutine RGYR ( \begin_inset LatexCommand ref reference "ite:Subroutine-RGYR" \end_inset )); the number of helical ( \begin_inset Formula $\beta$ \end_inset -sheet-like) residues 'nhel' ('nbet') and segments 'mhel' (`mbet') using subroutine HELIX ( \begin_inset LatexCommand ref reference "ite:subroutine-HELIX" \end_inset ); and the number of hydrogen bonds 'mhb' calculated through subroutine HBOND ( \begin_inset LatexCommand ref reference "ite:Subroutine-HBOND" \end_inset ). \end_layout \begin_layout Subsubsection Subroutine ANNEAL(nequi, nswp, nmes, tmax, tmin, lrand) \begin_inset LatexCommand index name "ANNEAL" \end_inset (file: anneal.f): \end_layout \begin_layout Standard (Usage: ' \emph on call anneal(nequi, nswp, nmes, tmax, tmin, lrand) \emph default ') \end_layout \begin_layout Standard This subroutine allows for optimization of protein conformations by means of simulated annealing. When the logical parameter 'lrand' is set to \shape italic .TRUE. \shape default then the program will start from a random configuration. Otherwise, the annealing process will start from the current configuration as stored in the array 'vlvr'. You have to set the initial and final temperatures 'tmax' and 'tmin', the number of Monte Carlo sweeps for simulation nswp and equilibration 'nequi', and the number 'nmes' of Monte Carlo sweeps between measurements. Four output files are created in the same directory where the executable is running: three of them contain the start (start.pdb), the final (final.pdb) and the global minimum configuration (global.pdb), respectively, the 4th file (time.d) contains the time series of the simulation. In the presented template we measure and store in this file the actual temperature 'temp', the total energy 'eol', the radius-of-gyration 'rgy' (as calculated by subroutine RGYR) and the Zimmerman code (as character variable 'zimm') of the present configuration as obtained through subroutine ZIMMER. \end_layout \begin_layout Subsubsection Subroutine MULCAN_PAR(nsweep, nup, temp, kmin, kmax, binWidth, l_iter) \begin_inset LatexCommand index name "MULCANPAR" \end_inset (file mulcan_mod_par.f): \end_layout \begin_layout Standard This subroutine is part of the module \emph on multicanonical \emph default . To use the function include the module with \emph on use multicanonical \emph default in your main file. The subroutine calculates multicanonical weight factors and has two working modes. One is for improving the already existing set of multicanonical parameters through further iterations, and another to calculate them anew. In the first case the logical variable 'l_iter' is set to the value \shape italic .TRUE. \shape default . The histogram from a preliminary simulation as well as the existing values of the parameters are read from the file 'mpar_full.d'. In the second case ( 'l_iter' = .FALSE.) all arrays are initialized to zero. The values of 'kmin' and 'kmax' together with and the size of the energy bin 'ebinWidth' determine the energy range . The total number of Monte Carlo sweeps is given by 'nsweep', 'nup' is the number of Monte Carlo sweeps between updates of the multicanonical parameters, and 'temp' is the temperature of the initial canonical simulation. Every 'nup' sweeps the histogram of the preliminary simulation as well as the existing values of the parameters are written to the file 'mpar_full.d'. The final multicanonical parameters are written into the file 'muca.d'. \end_layout \begin_layout Subsubsection Subroutine MULCAN_SIM(nequi, nsweeps, nmes, nsave, kmin, kmax, binWidth, restart) \begin_inset LatexCommand index name "MULCANSIM" \end_inset (file mulcan_sim.f): \end_layout \begin_layout Standard This subroutine performs a Monte Carlo simulation of a protein using multicanoni cal weights. Information necessary for a re-start are saved after every 'nsave' Monte Carlo moves into the file 'start.d'. You can restart a simulation by setting 'restart'=.TRUE. otherwise the simulation starts from a random configuration. Parameters 'binWidth', 'kmin', and 'kmax' set the energy-bin size and the lower and upper limits of the energy interval, within which the multicanonical parameters were defined (see MULCAN_PAR). These parameters are read from the file 'muca.d' and have to be the same as the ones used in MULCAN_PAR to generate the multicanonical weights. The total number of Monte Carlo sweeps is 'nsweep', the number of Monte Carlo sweeps needed for equilibration is 'nequi', and the interval between successive measurements is 'nmes'. In our example the output goes to 'time.d', in which the time series of the simulation is written. Measured quantities are the total potential energy 'eol', the number of contacts 'nhx' and the number of native contacts 'nhy'. The latter number counts the contacts that are the same in the present configuration and the target structure. Also measured is the Hamming distance 'dham' between the present and the target configuration. The latter three quantities are calculated by the subroutine CONTACTS. The contact map needs to be initialized before calling mulcan_sim., for example with the following code: \begin_inset ERT status collapsed \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard open(9,file='enkefa.ref') \end_layout \begin_layout Standard nresi=irsml2(1)-irsml1(1) + 1 \end_layout \begin_layout Standard ! nresi: Number of residues \end_layout \begin_layout Standard ! Read reference contact map \end_layout \begin_layout Standard nci = 0 \end_layout \begin_layout Standard do i=1,nresi \end_layout \begin_layout Standard read(9,*) (iref(i,j) , j=1,nresi) \end_layout \begin_layout Standard end do \end_layout \begin_layout Standard do i=1,nresi \end_layout \begin_layout Standard do j=nresi,i+3,-1 \end_layout \begin_layout Standard if(iref(i,j).eq.1) nci = nci + 1 \end_layout \begin_layout Standard end do \end_layout \begin_layout Standard end do \end_layout \begin_layout Standard write(*,*) 'Number of contacts in reference conformation:',nci \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Subsubsection Subroutine PARTEM_S \begin_inset LatexCommand index name "PARTEMS" \end_inset \end_layout \begin_layout Standard (file partem_s.f): \end_layout \begin_layout Standard This parameterless subroutine is designed to simulate a protein by means of the parallel tempering technique on a single processor. Before running a parallel tempering simulation using PARTEM_S one has to set the following variables in this subroutine: the number of replicas ('no') and the initial temperatures on each replica (array 'temp(no)'); the sum of the number of degrees of freedom over all replicas ('nvrmax'); the number of Monte Carlo sweeps at equilibration ('nequi') and the number of Monte Carlo sweeps between measurements ('nmes'). Except for re-starts, the logical variable 'lstart' should be set to .FALSE., in which case all dihedral angles are set to random values. At the end of the simulation all dihedral angles are stored in the file 'startconf.d', from which the program reads its start configuration, if 'lstart=.TRUE.'. In the simulation, measurements are made every 'nmes' Monte Carlo sweeps. In our template, we measure for each replica the radius-of gyration 'rgy' and the end-to-end distance 'ee' by means of subroutine RGYR and write both values together with the temperature 'temp0' , total energy 'energy', Monte Carlo time 'nsw', and the index of the replicas 'k1' to an output file 'time.d'. \end_layout \begin_layout Subsubsection PMAIN and subroutine PARTEM_P \begin_inset LatexCommand index name "PARTEMP" \end_inset (files pmain.f and partem_p.f) \end_layout \begin_layout Standard (Usage: ' \emph on call partem_p(no, nequi, nswp, nmes, newsta, switch) \emph default ') \end_layout \begin_layout Standard PMAIN and PARTEM_P are used for parallel tempering simulations on a parallel computer with 'no' nodes. PMAIN replaces MAIN in this case. Both PMAIN and PARTEM_P utilize the MPI message-passing routines. The number of nodes 'no' has to be set in the main program and equals the number of replicas in the parallel tempering method. In addition, one has to pass the following variables to PARTEM_P: the number of Monte Carlo sweeps for equilibration nequi and the total simulation nswp; and the number of Monte Carlo sweeps between measurements ('nmes'). Except for restarts the logical variable 'newsta' should be set to .FALSE.. If newsta is .TRUE., switch determines the starting configuration. If switch=1, all dihedral angles are set to the random values. If switch=-1, all dihedral angles are set to 180. If switch=0, the configuration remains unchanged. At the end of the simulation all dihedral angles are stored in the file 'par_R.in', from which the program reads its start configurations if 'newsta=.TRU E.'. Before running the program one has to prepare the file 'temp.d' with 'no' lines. Each line lists the index of replica 'j' and its temperature 'temp'. In our template we measure every 'nmes' Monte Carlo sweeps for each replica a number of quantities: the radius of gyration 'rgy', the number of residues 'nhel' and segments 'mhel' in an \begin_inset Formula $\alpha$ \end_inset -helix, the number of residues in a \begin_inset Formula $\beta$ \end_inset -sheet 'nbet' and the number of such segments 'mbet', the number of hydrogen bonds 'mhb', and the number of contacts 'nctot' and native contacts 'ncnat'. The latter number counts the contacts which are same in the present configurati on and in the target structure For this, one has to prepare in advance the file 'reference.d' with the contact matrix of the target structure. All measured quantities are written together with energy, temperature and the index of the corresponding replica in the file 'ts.d'. At the end of the simulation the average energy and specific heat for each temperature/replica and the acceptance rate of replica-exchange moves are written to standard output. \end_layout \begin_layout Subsection Useful Subroutines \end_layout \begin_layout Subsubsection Subroutine CONTACTS \begin_inset LatexCommand index name "CONTACTS" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic call contacts(ncn,nham2,dham) \shape default ): \end_layout \begin_layout Standard This subroutine calculates the matrix of the van-der-Waals contacts between C \begin_inset Formula $^{\alpha}$ \end_inset -atoms for a given configuration. 'ncn' is the total number of such contacts, 'nham2' the number of native contacts (i.e. such con- tacts which are the same in the present configuration and the reference configuration), and the Hamming distance 'dham' between the present and the target configuration. \end_layout \begin_layout Subsubsection Subroutine CNTENY \begin_inset LatexCommand index name "CNTENY" \end_inset : \end_layout \begin_layout Standard (Usage: \emph on call cnteny(nml) \emph default ) \end_layout \begin_layout Standard The subroutine ' \emph on cnteny \emph default ' displays a list of heavy atoms in molecule ' \emph on nml \emph default ' that have a contact energy of more than 2 kcal/mol. If one sets the parameter ' \emph on ieyel = 0 \emph default ' in ' \emph on cnteny \emph default ' the the contact energy is calculated with the van-der-Waals energy term, only. Setting ' \emph on ieyel = 1 \emph default ' additionally includes the electrostatic energy into the contact energy, and bad contacts are displayed if this energy exceeds a value of 10 kcal/mol for an atom. \end_layout \begin_layout Subsubsection Functions ENERGY \begin_inset LatexCommand index name "ENERGY" \end_inset , ENYSHE \begin_inset LatexCommand index name "ENYSHE" \end_inset , ENYFLX \begin_inset LatexCommand index name "ENYFLX" \end_inset , ENYSOL \begin_inset LatexCommand index name "ENYSOL" \end_inset , ESOLAN \begin_inset LatexCommand index name "ESOLAN" \end_inset , and ENYREG \begin_inset LatexCommand index name "ENYREG" \end_inset : \end_layout \begin_layout Standard (Usage: \shape italic eny = energy() \shape default ) \newline The real*8 function ` \shape italic energy \shape default ' is a wrapping function that calculates the energy of a given configuration. Depending on the choice of parameters ' \shape italic flex \shape default ', ' \shape italic sh2 \shape default ', ' \shape italic itysol \shape default ' and ' \shape italic ireg \shape default ' in ' \shape italic main \shape default ' (' \shape italic pmain \shape default ') the `vacuum' energy of the protein is calculated by the functions ' \shape italic enyshe \shape default ' (ECEPP/2 or ECEPP/3), or ' \shape italic enyflex \shape default ' (FLEX), and the solvation energy using the solvent accessible area method with ' \shape italic enysol \shape default ' (approximate, but fast estimation of the solvated area) or ' \shape italic esolan \shape default ' (analytical, but slower calculation of the area). ' \shape italic enyreg \shape default ' calculates the constraint energy term needed during regularization of PDB-structures. ' \shape italic energy \shape default ' returns the sum of these energy terms as selected by the parameters in the 'main'-module. \end_layout \begin_layout Subsubsection Subroutine GRADIENT \begin_inset LatexCommand index name "GRADIENT" \end_inset , OPESHE \begin_inset LatexCommand index name "OPESHE" \end_inset , OPEFLX \begin_inset LatexCommand index name "OPEFLX" \end_inset , OPESOL \begin_inset LatexCommand index name "OPESOL" \end_inset , and OPEREG \begin_inset LatexCommand index name "OPEREG" \end_inset : \end_layout \begin_layout Standard (Usage: \shape italic call gradient() \shape default ) \newline The wrapping subroutine ' \shape italic gradient \shape default ' calculates the energy and partial derivatives (with respect to internal degrees of freedom) for a given protein configuration. Depending on the choice of parameters ' \shape italic flex \shape default ', ' \shape italic sh2 \shape default ', ' \shape italic itysol \shape default ', and ' \shape italic ireg \shape default ', the subroutines ' \shape italic opeshe \shape default ' (ECEPP/2 or ECEPP/3), ' \shape italic opeflx \shape default ' (FLEX), ' \shape italic opesol \shape default ' (solvent accessible surface area) and/or ' \shape italic opereg \shape default ' (constraint energy term during regularization of PDB structures) are selected and contribute to the total energy and its gradient that is calculated. \end_layout \begin_layout Subsubsection \begin_inset LatexCommand label name "ite:Subroutine-HBOND" \end_inset Subroutine HBOND \begin_inset LatexCommand index name "HBOND" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic call hbond(nml,mhb,ipr) \shape default ): \newline This subroutine determines the number 'mhb' of the hydrogen bonds for the current conformation of the molecule and requires 'nml' as input - the index of the present molecule. The hydrogen bonds are printed out for ' \begin_inset Formula $ipr>0$ \end_inset '. \end_layout \begin_layout Subsubsection \begin_inset LatexCommand label name "ite:subroutine-HELIX" \end_inset subroutine HELIX \begin_inset LatexCommand index name "HELIX" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic helix(nhel,mhel,nbet,mbet) \shape default ): \newline HELIX determines whether a residue is part of an \begin_inset Formula $\alpha$ \end_inset -helix or \begin_inset Formula $\beta$ \end_inset -sheet according to its backbone dihedral angles ( \begin_inset Formula $\phi,\psi$ \end_inset ) (ranges defined in a PARAMETER statement) and returnes the number of helical residues 'nhel', the number of helical segments 'mhel', and the corresponding quantities for \begin_inset Formula $\beta$ \end_inset -strand elements 'nbet' and 'mbet'. \end_layout \begin_layout Subsubsection \begin_inset LatexCommand label name "ite:Subroutine-METROPOLIS" \end_inset Subroutine METROPOLIS \begin_inset LatexCommand index name "METROPOLIS" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic call metropolis(eol,acz,weight) \shape default ): \newline The subroutine implements one Metropolis move for every dihedral angle and requires as input a weight function ' \shape italic weight(eol) \shape default ' and the energy 'eol' of the present configuration. It returns in the variable 'eol' the energy of the configuration after 'nvr' (the number of non-fixed dihedral angles) updates. 'acz' accounts for the number of Metropolis moves were the proposed configurati on was accepted. \end_layout \begin_layout Subsubsection Subroutine OUTPDB \begin_inset LatexCommand index name "OUTPDB" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic call outpdb(nml,pdbFileName) \shape default ): \newline The subroutine writes in the output file the coordinates of the atoms of given molecule. The format is compatible with the standard data representation format of the Brookhaven Protein Data bank. This subroutine requires as input the index 'nml' of the given molecule and the file name of the output-file. If nml=0, the configuration for all molecules in the system is written to a multi-model pdb file. \end_layout \begin_layout Subsubsection Subroutine OUTVAR \begin_inset LatexCommand index name "OUTVAR" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic call outvar(nml, filename) \shape default ): \end_layout \begin_layout Standard This subroutine writes the dihedral angles of the current configuration of molecule ' \shape italic nml \shape default ' to either standard output ('filename = \begin_inset Quotes eld \end_inset \begin_inset Quotes eld \end_inset ') or a dedicated file, where filename is the name of the file). The dihedral angles are written out in a form that allows to use them as configuration input file for SMMP. If nml=0, the configuration for all molecules in the system is printed. \end_layout \begin_layout Subsubsection Function grnd() and subroutine sgrnd(seed): \end_layout \begin_layout Standard Our package includes an implementation of the twister-mersenne random generator. Initialize the random number generator using \emph on call sgrnd(seed) \emph default , where \emph on seed \emph default is an integer. Each call off \emph on grnd() \emph default returns a new random number between 0 and 1 taken from a uniform distribution. \end_layout \begin_layout Subsubsection \begin_inset LatexCommand label name "ite:Subroutine-RGYR" \end_inset Subroutine RGYR \end_layout \begin_layout Standard (Usage: \shape italic call rgyr(nml,rgy,ee) \shape default ): \newline This subroutine calculates the radius-of-gyration 'rgy' and the end-to-end distance 'ee' for the protein molecule 'nml'. Both quantities characterize the compactness of a protein conformation. \end_layout \begin_layout Subsubsection Function RMSDFUN \begin_inset LatexCommand index name "RMSDFUN" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic rmsd = rmsdfun(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl) \shape default ) \newline The real*8 function ' \shape italic rmsdfun \shape default ' calculates the root-mean square-deviation (rmsd) between atom postions in the current configuration of molecule ' \shape italic nml \shape default ' for a range of residues ' \shape italic [ir1,ir2] \shape default ' and a reference structure given by arrays of atom coordinates ' \shape italic xrf,yrf,zrf \shape default '. For each atom position in the current SMMP-structure the array ' \shape italic ixat \shape default ' provides the index of the equivalent atom in the reference structure, or yields a value of ' \shape italic 0 \shape default ' if there is no equivalent. By setting ' \shape italic isl=0 \shape default ' the rmsd is calculated over all non-hydrogen (heavy) atoms, for ' \shape italic isl = 1 \shape default ' for backbone atoms only, and for ' \shape italic isl = 2 \shape default ' the rmsd is calculated only between C \begin_inset Formula $^{\alpha}$ \end_inset -atoms. The routine ' \shape italic rmsdfun \shape default ' may only be called \series bold after \series default an initial call of the subroutine ' \shape italic rmsinit(nml,string) \shape default '. If \shape italic string = 'smmp' \shape default the reference configuration is a SMMP configuration, otherwise it is read in from a PDB-file which name given by ' \shape italic string \shape default '. \end_layout \begin_layout Subsubsection Subroutine ZIMMER \begin_inset LatexCommand index name "ZIMMER" \end_inset \end_layout \begin_layout Standard (Usage: \shape italic call zimmer(nresi) \shape default ): \newline This subroutine analyzes the backbone dihedrals of a given conformation and characterizes the state of each pair of ( \begin_inset Formula $\phi,\psi$ \end_inset ) by its Zimmerman code. The subroutine requires as input the number of residues 'nresi' and returns the character string 'zimm' with the Zimmerman code of the present configuratio n (through a common block). \end_layout \begin_layout Section How to cite SMMP \end_layout \begin_layout Standard Use of SMMP should be acknowledged by quoting the following reference: F. Eisenmenger, U.H.E. Hansmann, Sh. Hayryan, C.-K. Hu, \series bold [SMMP] A Modern Package for Simulation of Proteins \series default , \shape italic Comp. Phys. Comm. \shape default \series bold 138 \series default (2001) 192-212; \newline F. Eisenmenger, U.H.E. Hansmann, Sh. Hayryan, C.-K. Hu, \series bold An Enhanced Version of SMMP - Open-Source Software Package for Simulation of Proteins \series default , \shape italic Comp.Phys. Comm. 174 (2006) 422- \shape default 429; \newline \emph on J. H. Meinke, S. Mohanty, F. Eisenmenger and U.H.E. Hansmann, \series bold \emph default SMMP v. 3.0 --- Simulating proteins and protein interactions in Python and Fortran \series default , submitted to \emph on Comp.\InsetSpace ~ Phys.\InsetSpace ~ Comm. \newline \emph default Where appropriate, the program package's web-site should be quoted as: \series bold SMMP, \begin_inset LatexCommand url target "http://www.phy.mtu.edu/biophys/smmp.htm" \end_inset or \series default \begin_inset LatexCommand url target "http://www.smmp05.net" \end_inset . \end_layout \begin_layout Section Examples \begin_inset LatexCommand label name "sec:Examples" \end_inset \end_layout \begin_layout Standard All the examples can be found in the EXAMPLES subdirectory and should be run from there. All the serial examples can be built either from the top-level directory with \emph on make examples \emph default or from within the EXAMPLES directory with \emph on make \emph default . The MPI parallel tempering example is built from within the EXAMPLES directory with \emph on make parallel_tempering_p. \emph default The structure of all examples is very similar to Listing \begin_inset LatexCommand ref reference "listing:mainmc" \end_inset . We first initialize the force field we want to use, load the protein either from a sequence and a configuration file or from a PDB file, set the simulation parameters, and finally call the simulation subroutine. We use \shape italic \emph on \emph default grpn = 'nh2', grpc = 'cooh' \emph on as end groups in all of our examples. We also use a fixed value for the dielectric constant, so that \emph default epsd=.false. \emph on in all cases. \end_layout \begin_layout Standard The first example \begin_inset LatexCommand eqref reference "sub:Minimization-of-vacuum" \end_inset minimizes the energy of met-enkaphalin, a small endorphin with only 5 amino acids using a quasi-newton minimizer. The second examples \begin_inset LatexCommand eqref reference "sub:Structure-Determination-by" \end_inset tries to find the minimum vacuum energy structure of met-enkaphalin using simulated annealing. The third example \begin_inset LatexCommand eqref reference "sub:Calculation-of-Multicanonical" \end_inset shows how to calculate the estmated weights necessary for a multi-canonical Monte Carlo simulation. Met-enkaphalin again serves as sample system. In \begin_inset LatexCommand eqref reference "sub:Regularization-of-the" \end_inset we show how to regularize an experimental structure read from a PDB file. The force fields implemented use fixed bond lengths and valence angles. Without regularization, a structure read from a PDB file often has an unreasona bly high energy. Finally, we show how to perform parallel tempering simulations on multiple processors \begin_inset LatexCommand eqref reference "sub:Parallel-Tempering-simulation" \end_inset and on a single processor \begin_inset LatexCommand eqref reference "sub:Parallel-tempering-on" \end_inset . The sample output was generated on a system with two dual-core AMD Opteron 2216 processor running SUSE LINUX 10.1 (X86-64). The programs were compiled with Intels Fortran compiler with optimization set to -O2. \end_layout \begin_layout Subsection Minimization of vacuum energy from a given configuration \begin_inset LatexCommand label name "sub:Minimization-of-vacuum" \end_inset \end_layout \begin_layout Standard The MAIN module in minimization.f, has the following settings: \newline \shape italic \emph on ientyp = 0 \newline sh2 = .false. \newline itysol = 0 \newline iabin = 1 \shape default , \newline \emph default \newline which restrict energy calculations to evaluation of the ECEPP/3 energy (no protein-solvent interactions taken into account), choses \shape italic nh2 \shape default as N-terminal group and at the C-terminus the \shape italic cooh \shape default -group, and the input is read from a sequence/configuration file. The simulation task subroutine is called with \newline imin = 1 \newline maxit=10000 \newline eps=1.0d-9 \newline \shape italic call minim(imin,mait,eps) \shape default , \newline i.e. the simulation task is the local minimization of the vacuum energy of a configuration using the Quasi-Newton algorithm. Compile the program with \emph on make minimization \emph default , SMMP with this default settings will minimize the configuration given in \shape italic enkefa.var \shape default . Running minimization from the EXAMPLES directory will lead to the following output: \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard init_energy: itysol = 0 \end_layout \begin_layout Standard init_energy: esol_scaling = F \end_layout \begin_layout Standard init_molecule: Solvent: 0 \end_layout \begin_layout Standard File with sequence is enkefa.seq \end_layout \begin_layout Standard addend> Net charge of molecule # 1: 0.001 \end_layout \begin_layout Standard \end_layout \begin_layout Standard enkefa.var \end_layout \begin_layout Standard \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 1 Tyr : x1 set -172.590 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 1 Tyr : x2 set 78.710 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 1 Tyr : x6 set -165.880 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 1 Tyr : phi set -86.240 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 2 Gly : psi set 156.180 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 2 Gly : omg set -180.000 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 2 Gly : phi set -154.530 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 3 Gly : psi set 83.640 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 3 Gly : omg set 180.000 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 3 Gly : phi set 83.660 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 4 Phe : psi set -73.860 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 4 Phe : omg set -180.000 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 4 Phe : x1 set 58.790 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 4 Phe : x2 set 94.600 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 4 Phe : phi set -137.040 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : psi set 19.330 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : omg set -180.000 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : x1 set 52.760 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : x2 set 175.280 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : x3 set -179.830 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : x4 set -58.570 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : phi set -163.630 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : pst set 160.450 \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : omt set 180.000 \end_layout \begin_layout Standard \end_layout \begin_layout Standard \end_layout \begin_layout Standard Energy BEFORE minimization: \end_layout \begin_layout Standard \end_layout \begin_layout Standard Total: 0.10354E+04 \end_layout \begin_layout Standard Coulomb: 0.1991E+02 Lennard-Jones: 0.1142E+03 HB: 0.9008E+03 \end_layout \begin_layout Standard Variables: 0.4511E+00 Solvatation: 0.0000E+00 \end_layout \begin_layout Standard \end_layout \begin_layout Standard Step 1: energy 0.489976E+05 ( 0.867243E+13 ) \end_layout \begin_layout Standard Step 2: energy 0.108241E+03 ( 0.964364E+07 ) \end_layout \begin_layout Standard Step 3: energy -0.632755E+00 ( 0.498256E+04 ) \end_layout \begin_layout Standard Step 4: energy -0.739539E+01 ( 0.109979E+05 ) \end_layout \begin_layout Standard Step 5: energy 0.517957E+03 ( 0.678892E+09 ) \end_layout \begin_layout Standard Step 6: energy -0.767774E+01 ( 0.143880E+05 ) \end_layout \begin_layout Standard Step 7: energy -0.818570E+01 ( 0.100364E+05 ) \end_layout \begin_layout Standard Step 8: energy -0.938401E+01 ( 0.172508E+04 ) \end_layout \begin_layout Standard Step 9: energy -0.109799E+02 ( 0.105718E+04 ) \end_layout \begin_layout Standard Step 10: energy -0.116244E+02 ( 0.580821E+03 ) \end_layout \begin_layout Standard Step 11: energy -0.116969E+02 ( 0.145676E+04 ) \end_layout \begin_layout Standard Step 12: energy -0.118569E+02 ( 0.856882E+03 ) \end_layout \begin_layout Standard Step 13: energy -0.119576E+02 ( 0.521381E+03 ) \end_layout \begin_layout Standard Step 14: energy -0.120897E+02 ( 0.398602E+03 ) \end_layout \begin_layout Standard Step 15: energy -0.121476E+02 ( 0.222201E+03 ) \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard Step 70: energy -0.124285E+02 ( 0.194869E-04 ) \end_layout \begin_layout Standard Step 71: energy -0.124285E+02 ( 0.120524E-04 ) \end_layout \begin_layout Standard Step 72: energy -0.124285E+02 ( 0.641640E-05 ) \end_layout \begin_layout Standard Step 73: energy -0.124285E+02 ( 0.182141E-05 ) \end_layout \begin_layout Standard Step 74: energy -0.124285E+02 ( 0.219308E-05 ) \end_layout \begin_layout Standard Step 75: energy -0.124285E+02 ( 0.348349E-05 ) \end_layout \begin_layout Standard Step 76: energy -0.124285E+02 ( 0.166901E-05 ) \end_layout \begin_layout Standard Step 77: energy -0.124285E+02 ( 0.177102E-06 ) \end_layout \begin_layout Standard Step 78: energy -0.124285E+02 ( 0.663522E-08 ) \end_layout \begin_layout Standard Step 79: energy -0.124285E+02 ( 0.210693E-08 ) \end_layout \begin_layout Standard Step 80: energy -0.124285E+02 ( 0.370905E-08 ) \end_layout \begin_layout Standard Step 81: energy -0.124285E+02 ( 0.324300E-08 ) \end_layout \begin_layout Standard Step 82: energy -0.124285E+02 ( 0.777363E-09 ) \end_layout \begin_layout Standard Step 83: energy -0.124285E+02 ( 0.443358E-10 ) \end_layout \begin_layout Standard Step 84: energy -0.124285E+02 ( 0.322695E-11 ) \end_layout \begin_layout Standard Step 85: energy -0.124285E+02 ( 0.144369E-11 ) \end_layout \begin_layout Standard Step 86: energy -0.124285E+02 ( 0.291517E-11 ) \end_layout \begin_layout Standard Step 87: energy -0.124285E+02 ( 0.322695E-11 ) \end_layout \begin_layout Standard ---- CONVERGENCE ---- \end_layout \begin_layout Standard \end_layout \begin_layout Standard Final energies __________________________________________________ \end_layout \begin_layout Standard Total: -0.12429E+02 \end_layout \begin_layout Standard Coulomb: 0.2143E+02 Lennard-Jones: -0.2923E+02 HB: -0.6706E+01 \end_layout \begin_layout Standard Variables: 0.2084E+01 Solvatation: 0.0000E+00 \end_layout \begin_layout Standard \end_layout \begin_layout Standard Variables _________________ \end_layout \begin_layout Standard \end_layout \begin_layout Standard x1 1 -173.2 ( 0.6) \end_layout \begin_layout Standard x2 1 79.3 ( 0.6) \end_layout \begin_layout Standard x6 1 -166.3 ( 0.5) \end_layout \begin_layout Standard phi 1 -83.1 ( 3.2) \end_layout \begin_layout Standard psi 2 155.8 ( 0.4) \end_layout \begin_layout Standard omg 2 -177.1 ( 2.9) \end_layout \begin_layout Standard phi 2 -154.2 ( 0.3) \end_layout \begin_layout Standard psi 3 85.8 ( 2.2) \end_layout \begin_layout Standard omg 3 168.5 ( 11.5) \end_layout \begin_layout Standard phi 3 83.0 ( 0.7) \end_layout \begin_layout Standard psi 4 -75.0 ( 1.2) \end_layout \begin_layout Standard omg 4 -170.0 ( 10.0) \end_layout \begin_layout Standard x1 4 58.9 ( 0.1) \end_layout \begin_layout Standard x2 4 94.5 ( 0.1) \end_layout \begin_layout Standard phi 4 -136.8 ( 0.2) \end_layout \begin_layout Standard psi 5 19.1 ( 0.2) \end_layout \begin_layout Standard omg 5 -174.1 ( 5.9) \end_layout \begin_layout Standard x1 5 52.9 ( 0.1) \end_layout \begin_layout Standard x2 5 175.3 ( 0.0) \end_layout \begin_layout Standard x3 5 -179.9 ( 0.0) \end_layout \begin_layout Standard x4 5 -58.6 ( 0.0) \end_layout \begin_layout Standard phi 5 -163.4 ( 0.2) \end_layout \begin_layout Standard pst 5 160.8 ( 0.3) \end_layout \begin_layout Standard omt 5 -179.8 ( 0.2) \end_layout \begin_layout Standard \end_layout \begin_layout Standard Gradient______________________________________________________________ \end_layout \begin_layout Standard 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 \end_layout \begin_layout Standard 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 \end_layout \begin_layout Standard 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Subsection Structure Determination by Simulated Annealing \begin_inset LatexCommand label name "sub:Structure-Determination-by" \end_inset \end_layout \begin_layout Standard For the following simulated annealing run for Met-Enkephalin defined through its sequence file \shape italic enkefa.seq \shape default , another configuration file \shape italic enkefa.ann \shape default is used. In this configuration file the dihedral angles \shape italic omg \shape default and \shape italic omt \shape default are fixed at 180.0 degrees using \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard om*: 180.0 & \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset . Simulated annealing performs only poorly when these angles are free to change. In this example we choose ECEPP/2 ( setting \shape italic sh2 = .true. \shape default and no protein-solvation interaction. Hence, we have the following settings for parameters in \shape italic annealing.f \shape default : \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard ientyp = 0 \end_layout \begin_layout Standard sh2 = .true. \end_layout \begin_layout Standard itysol = 0 \end_layout \begin_layout Standard iabin = 1 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset . Our simulated annealing run consists of 100 Monte Carlo sweeps at the highest temperature (1000 K) followed by 100,000 sweeps during which the temperature is gradually decreased to a final temperature of 100 K. The progress of the simulation is monitored by writing data into the file ` \shape italic time.d \shape default ' every 1000 sweeps. The simulation starts from a random configuration. \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard lrand=.true. \end_layout \begin_layout Standard nequi=100, \end_layout \begin_layout Standard nswp=100000 \end_layout \begin_layout Standard nmes=1000 \end_layout \begin_layout Standard tmax=1000.0 \end_layout \begin_layout Standard tmin=100.0 \end_layout \begin_layout Standard \end_layout \begin_layout Standard call anneal(nequi, nsweeps, nmes, tmax, tmin, lrand) \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard After compiling the example with \emph on make annealing \emph default in the sub-directory \shape italic EXAMPLES, \shape default \shape italic \emph on we run the program by calling \shape default \emph default \shape italic annealing. \shape default The output will contain the energy of the start configuration, energy after equilibration, the average acceptance rate of moves in the simulated annealing run, the final configuration (and its energy) and the configuration with the lowest energy found in this run: \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard init_energy: itysol = 0 \end_layout \begin_layout Standard init_energy: esol_scaling = F \end_layout \begin_layout Standard init_molecule: Solvent: 0 \end_layout \begin_layout Standard File with SEQUENCE: enkefa.seq \end_layout \begin_layout Standard enkefa.ann \end_layout \begin_layout Standard \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 2 Gly : omg set 180.000 Fixed \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 3 Gly : omg set 180.000 Fixed \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 4 Phe : omg set 180.000 Fixed \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : omg set 180.000 Fixed \end_layout \begin_layout Standard redvar> Met-Enkephalin: residue 5 Met : omt set 180.000 Fixed \end_layout \begin_layout Standard redvar> Molecule Met-Enkephalin: 19 variable(s) remain unchanged \end_layout \begin_layout Standard \end_layout \begin_layout Standard energy of start configuration: 0.75921E+09 \end_layout \begin_layout Standard \end_layout \begin_layout Standard Energy after equilibration: 11.0835755006373 \end_layout \begin_layout Standard acceptance rate: 0.181430817270775 \end_layout \begin_layout Standard \end_layout \begin_layout Standard last energy -8.03108763810183 \end_layout \begin_layout Standard 1 : 1 : x1 : -168.428 \end_layout \begin_layout Standard 1 : 1 : x2 : 89.562 \end_layout \begin_layout Standard 1 : 1 : x6 : -2.861 \end_layout \begin_layout Standard 1 : 1 : phi : -112.654 \end_layout \begin_layout Standard 1 : 2 : psi : 156.420 \end_layout \begin_layout Standard 1 : 2 : omg : -180.000 & \end_layout \begin_layout Standard 1 : 2 : phi : -174.567 \end_layout \begin_layout Standard 1 : 3 : psi : 146.718 \end_layout \begin_layout Standard 1 : 3 : omg : -180.000 & \end_layout \begin_layout Standard 1 : 3 : phi : 68.148 \end_layout \begin_layout Standard 1 : 4 : psi : -96.563 \end_layout \begin_layout Standard 1 : 4 : omg : -180.000 & \end_layout \begin_layout Standard 1 : 4 : x1 : 59.864 \end_layout \begin_layout Standard 1 : 4 : x2 : -83.907 \end_layout \begin_layout Standard 1 : 4 : phi : -141.914 \end_layout \begin_layout Standard 1 : 5 : psi : 20.747 \end_layout \begin_layout Standard 1 : 5 : omg : -180.000 & \end_layout \begin_layout Standard 1 : 5 : x1 : 61.133 \end_layout \begin_layout Standard 1 : 5 : x2 : -174.118 \end_layout \begin_layout Standard 1 : 5 : x3 : 176.917 \end_layout \begin_layout Standard 1 : 5 : x4 : -178.096 \end_layout \begin_layout Standard 1 : 5 : phi : -162.246 \end_layout \begin_layout Standard 1 : 5 : pst : 161.616 \end_layout \begin_layout Standard 1 : 5 : omt : -180.000 & \end_layout \begin_layout Standard \end_layout \begin_layout Standard lowest energy ever found: 63297 -9.40635361798398 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset The random start configuration, final configuration and lowest-energy configura tion are also written in PDB-format into the files \shape italic start.pdb \shape default , \shape italic final.pdb \shape default and \shape italic global.pdb \shape default allowing for easy visualization of these configurations with programs graphical ly displaying PDB-structures. The progress of the simulated annealing run is monitored through the time series in file \shape italic time.d \shape default which list the Monte Carlo sweep, temperature, ECEPP/2 energy, radius-of-gyrati on, the partial energy terms \begin_inset Formula $E_{\mbox{hb}}$ \end_inset , \begin_inset Formula $E_{\mbox{vdW}}$ \end_inset , \begin_inset Formula $E_{\mbox{C}}$ \end_inset , \begin_inset Formula $E_{\mbox{vr}}$ \end_inset and the Zimmerman-coding of the present configuration. All the terms related to the solvent have been replaced by a single dot per column in the sample below. \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard # $Id: anneal.f 192 2007-01-02 14:06:25Z meinke $ \end_layout \begin_layout Standard # nsw, temp, eol, eysl, eyslh, eyslp, asa, rgy, \end_layout \begin_layout Standard # rgyh, rgyp, eyhb, eyvw, eyel, eyvr, zimm \end_layout \begin_layout Standard 0 1000.000 8.045 .... 4.542 .. -3.992 -14.929 23.308 3.658 ghcBD \end_layout \begin_layout Standard 1000 977.237 16.397 .... 5.835 .. -2.029 -11.821 24.061 6.186 BDeCG \end_layout \begin_layout Standard 2000 954.993 16.147 .... 5.302 .. -1.786 -11.072 25.230 3.775 AdEDA \end_layout \begin_layout Standard 3000 933.254 26.320 .... 6.434 .. -1.324 -4.642 25.514 6.771 EecFC \end_layout \begin_layout Standard 4000 912.011 14.884 .... 5.145 .. -2.474 -13.559 25.449 5.469 EHDDB \end_layout \begin_layout Standard 5000 891.251 21.650 .... 5.882 .. -1.179 -8.380 27.125 4.083 EHfEB \end_layout \begin_layout Standard 6000 870.964 13.854 .... 5.806 .. -2.271 -12.703 25.621 3.207 EcaEh \end_layout \begin_layout Standard 7000 851.138 19.786 .... 5.439 .. -1.952 -9.254 26.428 4.565 EHaAE \end_layout \begin_layout Standard 8000 831.764 14.299 .... 5.986 .. -1.596 -16.538 24.699 7.735 EECEF \end_layout \begin_layout Standard 9000 812.831 7.773 .... 4.915 .. -2.819 -16.392 24.850 2.134 EECBB \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard 91000 123.027 -8.381 .... 4.414 .. -5.162 -26.932 22.230 1.483 FEcDE \end_layout \begin_layout Standard 92000 120.226 -6.302 .... 4.439 .. -4.825 -25.727 22.824 1.426 EEcDE \end_layout \begin_layout Standard 93000 117.490 -6.655 .... 4.413 .. -4.355 -25.081 21.937 0.845 FEcDE \end_layout \begin_layout Standard 94000 114.815 -8.190 .... 4.523 .. -6.848 -24.025 21.598 1.086 FEcDE \end_layout \begin_layout Standard 95000 112.202 -5.965 .... 4.434 .. -5.052 -24.021 21.995 1.114 EEcDE \end_layout \begin_layout Standard 96000 109.648 -7.081 .... 4.444 .. -4.717 -25.137 22.669 0.104 EEcDE \end_layout \begin_layout Standard 97000 107.152 -8.042 .... 4.561 .. -4.969 -25.544 21.850 0.621 FEcDE \end_layout \begin_layout Standard 98000 104.713 -7.342 .... 4.439 .. -4.737 -26.152 22.463 1.083 FEcDE \end_layout \begin_layout Standard 99000 102.329 -6.561 .... 4.438 .. -3.975 -25.213 21.963 0.664 FEcDE \end_layout \begin_layout Standard 100000 100.000 -8.031 .... 4.433 .. -5.231 -25.792 22.659 0.333 EEcDE \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Subsection Calculation of Multicanonical Parameters \begin_inset LatexCommand label name "sub:Calculation-of-Multicanonical" \end_inset \end_layout \begin_layout Standard In the following example, multicanonical parameters are calculated for met-enkep halin (with fixed omega angles) in gas-phase simulations relying on the ECEPP/3 force-field. \shape italic \emph on The functions \shape default \emph default \shape italic mulcan_par \shape default \shape italic \emph on and \shape default \emph default \shape italic mulcan_sim \shape default are part of the module \emph on multicanonical. \emph default To use the module, include \emph on use multicanonical \emph default in \emph on main \emph default . For this example, we set the following parameters in ' \shape italic main \shape default ': \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard ientyp = 0 \end_layout \begin_layout Standard sh2 = .false. \end_layout \begin_layout Standard itysol = 0 \end_layout \begin_layout Standard iabin = 1 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard The simulation has a total of 100,000 Monte Carlo sweeps and the multicanonical parameters are updated every 2,000 sweeps. The parameters are calculated from scratch, and we attempt to obtain a histogram in energy in the interval [-12,20] Kcal/mol with energy bin size 1 Kcal/mol. We expect that at temperature 1000 K we are clearly in the high-energy phase where successive configurations have little correlation. This leads to the following set of arguments for calling \emph on mulcan \emph default _ \emph on par \emph default : \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard l_iter=.false. \end_layout \begin_layout Standard kmin=-12 \end_layout \begin_layout Standard kmax=20 \end_layout \begin_layout Standard ebin=1.0d0 \end_layout \begin_layout Standard nsweep=100000 \end_layout \begin_layout Standard nup=5000 \end_layout \begin_layout Standard temp=1000.0 \end_layout \begin_layout Standard call mulcan_par(nsweep, nup, temp, kmin, kmax, l_iter) \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard This example is available in the file \emph on multicanonical.f \emph default in the EXAMPLES directory. After compiling it with \emph on make \emph default \emph on multicanonical \emph default and running the program using \emph on multicanonical \emph default the calculated multicanonical parameters are written into the file ' \shape italic muca.d \shape default '. For further iteration of these parameters additional quantities such as the accumulated histogram are written into a file ' \shape italic mpar_full.d \shape default ' (out of which ' \shape italic mulcan_par \shape default ' reads them if ' \shape italic l_iter=.true. \shape default ' is set). We show here only the file ' \shape italic muca.d \shape default ' which lists the index of the energy bin, the \begin_inset Formula $\beta(E)$ \end_inset and the \begin_inset Formula $\alpha(E)$ \end_inset parameters (see \begin_inset LatexCommand citep key "Eisenmenger2001" \end_inset ) for more detailed explanations and references on multicanonical simulations): \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard -12 11.8671041270715 93.0982188015155 \end_layout \begin_layout Standard -11 6.74942546567488 34.2449141954539 \end_layout \begin_layout Standard -10 4.76239638525333 13.3811088510276 \end_layout \begin_layout Standard -9 2.62566106761319 -6.91787666655369 \end_layout \begin_layout Standard -8 2.76891298156071 -5.70023539799973 \end_layout \begin_layout Standard -7 2.98636199792932 -4.06936777523520 \end_layout \begin_layout Standard -6 2.62458524605883 -6.42091666239336 \end_layout \begin_layout Standard -5 2.68953283146133 -6.06370494267963 \end_layout \begin_layout Standard -4 1.99405651648455 -9.19334836007513 \end_layout \begin_layout Standard -3 2.37207855308226 -7.87027123198313 \end_layout \begin_layout Standard -2 1.80430850778520 -9.28969634522579 \end_layout \begin_layout Standard -1 1.91571656031594 -9.12258426642968 \end_layout \begin_layout Standard 0 1.72752866077015 -9.21667821620257 \end_layout \begin_layout Standard 1 1.66881722525666 -9.18732249844583 \end_layout \begin_layout Standard 2 1.50004119691055 -8.93415845592665 \end_layout \begin_layout Standard 3 1.39142908670518 -8.66262818041324 \end_layout \begin_layout Standard 4 1.16605491925085 -7.87381859432308 \end_layout \begin_layout Standard 5 1.88002935687926 -11.0867035636509 \end_layout \begin_layout Standard 6 0.784923700785834 -5.06362245513708 \end_layout \begin_layout Standard 7 1.15928001690741 -7.49693850992732 \end_layout \begin_layout Standard 8 1.53630835441332 -10.3246510412217 \end_layout \begin_layout Standard 9 0.804876443148135 -4.10747979546758 \end_layout \begin_layout Standard 10 1.20439627142663 -7.90291816411324 \end_layout \begin_layout Standard 11 0.843296201282429 -4.11136742759918 \end_layout \begin_layout Standard 12 1.03327436670481 -6.29611632995659 \end_layout \begin_layout Standard 13 0.847836723791197 -3.97814579353639 \end_layout \begin_layout Standard 14 0.794466332096526 -3.25764550565833 \end_layout \begin_layout Standard 15 0.665347460625295 -1.38542186932548 \end_layout \begin_layout Standard 16 0.897108190646473 -4.97771318465375 \end_layout \begin_layout Standard 17 0.527634627532446 1.11860060672770 \end_layout \begin_layout Standard 18 0.591554662202600 0.000000000000000E+000 \end_layout \begin_layout Standard 19 0.591554662202600 0.000000000000000E+000 \end_layout \begin_layout Standard 20 0.591554662202600 0.000000000000000E+000 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard A test simulation of 100,000 Monte Carlo sweeps (using subroutine ' \shape italic mulcan_sim \shape default ' with these weights led to the following histogram that is flat over the studied energy range (as is characteristic for the multicanonical ensemble) (note that the last energy bin contains all entries for energies larger/equal than 20 Kcal/mol): \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard -9 105 \end_layout \begin_layout Standard -8 828 \end_layout \begin_layout Standard -7 1945 \end_layout \begin_layout Standard -6 3108 \end_layout \begin_layout Standard -5 3247 \end_layout \begin_layout Standard -4 3054 \end_layout \begin_layout Standard -3 2844 \end_layout \begin_layout Standard -2 2531 \end_layout \begin_layout Standard -1 2397 \end_layout \begin_layout Standard 0 2082 \end_layout \begin_layout Standard 1 2179 \end_layout \begin_layout Standard 2 2294 \end_layout \begin_layout Standard 3 2579 \end_layout \begin_layout Standard 4 3193 \end_layout \begin_layout Standard 5 2855 \end_layout \begin_layout Standard 6 2927 \end_layout \begin_layout Standard 7 3994 \end_layout \begin_layout Standard 8 3362 \end_layout \begin_layout Standard 9 3047 \end_layout \begin_layout Standard 10 3289 \end_layout \begin_layout Standard 11 3199 \end_layout \begin_layout Standard 12 3412 \end_layout \begin_layout Standard 13 3034 \end_layout \begin_layout Standard 14 3092 \end_layout \begin_layout Standard 15 3258 \end_layout \begin_layout Standard 16 3079 \end_layout \begin_layout Standard 17 2957 \end_layout \begin_layout Standard 18 3157 \end_layout \begin_layout Standard 19 3055 \end_layout \begin_layout Standard 20 19897 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Subsection Regularization of the B domain of protein A \begin_inset LatexCommand label name "sub:Regularization-of-the" \end_inset \end_layout \begin_layout Standard Frequently, one wants to compare the results of a simulation with experimental data. When we compare properties of calculated configurations with those of structure s from the Protein Data Bank (PDB), we have to remember that the actual bond lengths and angles in the PDB-structure will differ from the fixed values that are assumed with the ECEPP, FLEX, and Lund potentials. Forcing the molecule into the standard bonding geometry model may lead to un-physically high energies. The process of finding an optimal structure within the standard geometry model starting from a PDB-structure is called regularization. \end_layout \begin_layout Standard In our example, we regularize of the coordinates for residues 10-55 of the B domain of \shape italic Staphylococcus Aureus \shape default Protein A, which were taken from entry \shape italic 1BDD \shape default of the Protein Data Bank (see \begin_inset LatexCommand url target "http://www.rcsb.org" \end_inset ), as provided with this package in file \shape italic 1bdd.pdb \shape default . Setting the variable ' \shape italic iabin = 0 \shape default ' in \shape italic main.f \shape default makes \shape italic init_molecule \shape default call subroutine \shape italic pdbread \shape default that reads the amino acid sequence and the atomic coordinates from the PDB- structure and stores this data in arrays declared in ' \shape italic INCP.H \shape default '. Then \shape italic pdbvars \shape default is automatically called to measure all dihedral angles. In order to regularize the PDB-structure one has to \shape italic call regul(nml, iter, nsteps, \begin_inset Formula $\varepsilon$ \end_inset ) \shape default as the task subroutine in ' \shape italic main \shape default ' . Set the following parameters in ' \shape italic main \shape default ', accordingly: \newline \newline \shape italic ientyp=0 \newline sh2 = .false. \newline itysol = 0 \newline iabin = 0 \shape default \newline \newline Note that ' \shape italic regul \shape default ', as any other task subroutine, has to be called AFTER ' \shape italic init_molecule \shape default ' . The example is available in the EXAMPLES directory. You can compile it by calling \emph on make regularization. \emph default Run the program by calling \emph on regularization \emph default . We use 10 iterations with up to 15000 sweeps each and \begin_inset Formula $\varepsilon=10^{-7}$ \end_inset . A value for \begin_inset Formula $\varepsilon$ \end_inset that is too large may lead to numerical instabilities. \end_layout \begin_layout Standard First ' \shape italic regul \shape default ' obtains its input from the arrays declared in ' \shape italic INCP.H \shape default '. In a first step a naive representation of the molecule is build in the ECEPP geometry from the stored dihedral angles that were calculated from the PDB coordinates. This simple model serves as starting configuration for the regularization. In our example, its root-mean-square deviation (over all heave atoms) compared to the PDB-structure initially is 2.8 \begin_inset ERT status collapsed \begin_layout Standard \backslash AA \end_layout \end_inset and its ECEPP-energy \begin_inset Formula $\approx10^{5}$ \end_inset kcal/mol, mainly due to excessive van-der-Waals repulsions. The regularization starts by first minimizing a term that measures the sum of squared distances between heavy atoms in the SMMP-structure and the given PDB-structure. This leads to a rmsd of 0.14 \begin_inset ERT status collapsed \begin_layout Standard \backslash AA \end_layout \end_inset . After this initial step a list of bad contacts (vdW-energy of more than 2 kcal/mol) in the SMMP-structure is printed out. In a second step the physical energy is minimized, allowing only the free hydrogens to move, and the bad contacts and rmsd are displayed again. This step reduces the ECEPP-energy to \begin_inset Formula $\approx10^{3}$ \end_inset kcal/mol. ' \shape italic Regul \shape default ' aims at further reducing this energy while at the same time keeping the rmsd as small as possible. This is done by minimizing a composite energy from the weighted sum of physical ECEPP-energy and constraint energy (the quadratic distance measure) over 10 iterations. In each iteration the weight of the constraint energy is successively lowered (from 1 to 0) and the weight of the ECEPP energy raised (from 0 to 1). At the end, the dihedral angles of the final structure are printed out together with its rmsd and a list of (eventually) remaining bad contacts. \end_layout \begin_layout Standard The following output is obtained from this last iteration with a weight of 1 for the ECEPP-energy and a weight of 0 for the constraint energy, that leads to a SMMP-configuration with a total ECEPP-energy of \begin_inset Formula $-429.9$ \end_inset kcal/mol and a rmsd of 1.89 \begin_inset ERT status inlined \begin_layout Standard \backslash AA \end_layout \end_inset , compared to the PDB-structure: \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard ================ Minimization #10 Wt(energy) = 0.100E+01 Wt(regul.) = 0.000E+00 \end_layout \begin_layout Standard \end_layout \begin_layout Standard \end_layout \begin_layout Standard Energy BEFORE minimization: \end_layout \begin_layout Standard \end_layout \begin_layout Standard Total: -0.39831E+03 \end_layout \begin_layout Standard Coulomb: -0.4685E+02 Lennard-Jones: -0.3306E+03 HB: -0.6915E+02 \end_layout \begin_layout Standard Variables: 0.4831E+02 Solvatation: 0.0000E+00 Regularization: 0.3315E+03 \end_layout \begin_layout Standard \end_layout \begin_layout Standard Step 1: energy 0.259565E+13 ( 0.296722E+30, 0.363456E+12 ) \end_layout \begin_layout Standard Step 2: energy 0.481884E+09 ( 0.572102E+22, 0.644342E+12 ) \end_layout \begin_layout Standard Step 3: energy 0.444100E+10 ( 0.599018E+24, 0.617413E+12 ) \end_layout \begin_layout Standard Step 4: energy 0.491253E+10 ( 0.220900E+25, 0.533540E+12 ) \end_layout \begin_layout Standard Step 5: energy -0.355950E+03 ( 0.341264E+06, 0.495213E+10 ) \end_layout \begin_layout Standard Step 6: energy -0.395467E+03 ( 0.762827E+06, 0.726715E+08 ) \end_layout \begin_layout Standard Step 7: energy -0.398489E+03 ( 0.198935E+05, 0.137122E+08 ) \end_layout \begin_layout Standard Step 8: energy -0.398761E+03 ( 0.409017E+05, 0.187175E+08 ) \end_layout \begin_layout Standard Step 9: energy -0.399231E+03 ( 0.197834E+05, 0.329551E+08 ) \end_layout \begin_layout Standard Step 10: energy -0.399787E+03 ( 0.143934E+06, 0.121392E+09 ) \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard Step 415: energy -0.429868E+03 ( 0.135007E-03, 0.523386E+11 ) \end_layout \begin_layout Standard Step 416: energy -0.429868E+03 ( 0.123681E-04, 0.523385E+11 ) \end_layout \begin_layout Standard Step 417: energy -0.429868E+03 ( 0.254870E-05, 0.523385E+11 ) \end_layout \begin_layout Standard Step 418: energy -0.429868E+03 ( 0.193334E-05, 0.523385E+11 ) \end_layout \begin_layout Standard Step 419: energy -0.429868E+03 ( 0.140748E-05, 0.523385E+11 ) \end_layout \begin_layout Standard Step 420: energy -0.429868E+03 ( 0.740869E-06, 0.523384E+11 ) \end_layout \begin_layout Standard Step 421: energy -0.429868E+03 ( 0.489693E-06, 0.523384E+11 ) \end_layout \begin_layout Standard Step 422: energy -0.429868E+03 ( 0.282056E-06, 0.523384E+11 ) \end_layout \begin_layout Standard Step 423: energy -0.429868E+03 ( 0.104563E-06, 0.523384E+11 ) \end_layout \begin_layout Standard Step 424: energy -0.429868E+03 ( 0.519142E-07, 0.523385E+11 ) \end_layout \begin_layout Standard Step 425: energy -0.429868E+03 ( 0.104563E-06, 0.523384E+11 ) \end_layout \begin_layout Standard ---- CONVERGENCE ---- \end_layout \begin_layout Standard \end_layout \begin_layout Standard Final energies __________________________________________________ \end_layout \begin_layout Standard \end_layout \begin_layout Standard Total: -0.42987E+03 \end_layout \begin_layout Standard Coulomb: -0.5431E+02 Lennard-Jones: -0.3510E+03 HB: -0.7268E+02 \end_layout \begin_layout Standard Variables: 0.4809E+02 Solvatation: 0.0000E+00 Regularization: 0.2241E+05 \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard \end_layout \begin_layout Standard RMSD = 1.89099626258629 \end_layout \begin_layout Standard \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Subsection Parallel Tempering simulation of four A \begin_inset Formula $\beta_{16-22}$ \end_inset peptides \begin_inset LatexCommand label name "sub:Parallel-Tempering-simulation" \end_inset \end_layout \begin_layout Standard In the following example, a parallel tempering simulation of the fragment 16--22 of Alzheimer's \begin_inset Formula $\beta$ \end_inset -amyloid in solvent. It relies on the ECEPP/3 force-field and all dihedral angles are free. References for the algorithm can be found in \begin_inset LatexCommand cite key "Eisenmenger2001" \end_inset . The simulation is performed on a parallel computer on \begin_inset Formula $n\times8$ \end_inset processors. The example is available in the EXAMPLES directory as \emph on parallel_tempering_p \emph default . Commpile it using \emph on make parallel_tempering \emph default . Set the following parameters \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard num_replica = 8 \end_layout \begin_layout Standard ientyp = 0 \end_layout \begin_layout Standard sh2=.false. \end_layout \begin_layout Standard itysol=1 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard After initialization, the program calls \emph on init_molecule \emph default once for every molecule. Only the last call contains a reference to a variable file. It contains the global coordinates for each molecule \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard @ 1 : 100, 100, 100, 0, 0, 0 \end_layout \begin_layout Standard @ 2 : -100, -100, -100, 0, 0, 0 \end_layout \begin_layout Standard @ 3 : 100, -100, 100, 0, 0, 0 \end_layout \begin_layout Standard @ 4 : 100, -100, -100, 0, 0, 0 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard The subroutine reads the distribution of temperatures (in our example 'temperatu res_abeta') that has to be provided. \end_layout \begin_layout Standard The simulation starts from a stretched configuration with 100 sweeps for equilibration. The molecule is then studied over 10,000 sweeps on each node, with measurements taken every 10 sweep. \end_layout \begin_layout Standard Hence the following arguments have to be passed \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard nequi=100 \end_layout \begin_layout Standard nswp=10000 \end_layout \begin_layout Standard nmes=10 \end_layout \begin_layout Standard nsave=1000 \end_layout \begin_layout Standard switch=0 \end_layout \begin_layout Standard newsta=.true. \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset by calling \emph on partem_p(num_replica, nequi, nswp, nmes, nsave, newsta, switch, rep_id, partem_comm). \end_layout \begin_layout Standard After compiling with \emph on make parallel_tempering \emph default the program is ready to run. Note that you must first execute \emph on make parallel \emph default in the top-level directory. We do not provide a script to launch the program since launching it depends on the specific MPI installation. The data are written into the file ' \shape italic ts.d \shape default '. In our case, the data are the sweep, the temperature index, the index of the replica this temperature is currently associated, the inverse temperature, and energy, radius-of gyration, number of residues that are part of a helix or sheet, number of hydrogen bond, total number of contacts and the number of contacts that appear also in the reference configuration: \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard 10 1 1 0.503086435280445 \end_layout \begin_layout Standard 18.5259057594418 -101.020909189261 1.668390368424223E-004 \end_layout \begin_layout Standard 153.870626474863 3 4 0 \end_layout \begin_layout Standard 10 2 2 0.838477392134076 \end_layout \begin_layout Standard -21.8597428563679 -112.839132565408 4.398819716124160E-005 \end_layout \begin_layout Standard 164.067677362049 0 12 0 \end_layout \begin_layout Standard 10 3 3 1.25771608820111 \end_layout \begin_layout Standard -63.6878966409724 -115.697268416431 7.324283587028595E-005 \end_layout \begin_layout Standard 162.383801785309 2 9 0 \end_layout \begin_layout Standard 10 4 4 1.43738981508699 \end_layout \begin_layout Standard -82.2329362395502 -118.850435462771 -2.589501558808733E-005 \end_layout \begin_layout Standard 155.367430914339 3 8 0 \end_layout \begin_layout Standard 10 5 5 1.67695478426815 \end_layout \begin_layout Standard -89.8268073782128 -115.511880821486 -1.910806982598894E-004 \end_layout \begin_layout Standard 156.197401324212 1 13 0 \end_layout \begin_layout Standard 10 6 6 1.82940521920162 \end_layout \begin_layout Standard -102.200054103571 -119.664523580541 -2.902461830424936E-005 \end_layout \begin_layout Standard 158.377278261187 1 7 0 \end_layout \begin_layout Standard 10 7 7 1.92017723389483 \end_layout \begin_layout Standard -92.5898004594444 -111.866857243732 8.412738265616228E-005 \end_layout \begin_layout Standard 156.691206902137 2 13 0 \end_layout \begin_layout Standard 10 8 8 2.01234574112178 \end_layout \begin_layout Standard -91.7374330640579 -118.481236403476 -1.450348059239640E-004 \end_layout \begin_layout Standard 160.611273802438 0 9 0 \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard ... \end_layout \begin_layout Standard 10000 1 2 0.503086435280445 \end_layout \begin_layout Standard 25.7121801637793 -112.113204586669 -7.981375910951304E-005 \end_layout \begin_layout Standard 171.130731978089 0 5 0 \end_layout \begin_layout Standard 10000 2 4 0.838477392134076 \end_layout \begin_layout Standard -83.1513580270155 -119.884471326212 -3.185178277838820E-005 \end_layout \begin_layout Standard 179.222034577900 4 5 0 \end_layout \begin_layout Standard 10000 3 3 1.25771608820111 \end_layout \begin_layout Standard -88.7115772801244 -109.660390359632 2.816222447889407E-005 \end_layout \begin_layout Standard 169.412616033656 6 2 0 \end_layout \begin_layout Standard 10000 4 1 1.43738981508699 \end_layout \begin_layout Standard -117.230892149812 -118.180458392844 -1.362571932424056E-004 \end_layout \begin_layout Standard 155.662046672612 9 3 0 \end_layout \begin_layout Standard 10000 5 5 1.67695478426815 \end_layout \begin_layout Standard -117.452066538358 -114.923193931984 9.292455852280616E-005 \end_layout \begin_layout Standard 165.838174859018 8 4 0 \end_layout \begin_layout Standard 10000 6 8 1.82940521920162 \end_layout \begin_layout Standard -114.553279087787 -116.322737335011 1.354655604376923E-005 \end_layout \begin_layout Standard 161.094515087733 4 4 0 \end_layout \begin_layout Standard 10000 7 7 1.92017723389483 \end_layout \begin_layout Standard -128.613090309274 -113.083234658370 6.618583172642824E-005 \end_layout \begin_layout Standard 169.226190038938 10 0 0 \end_layout \begin_layout Standard 10000 8 6 2.01234574112178 \end_layout \begin_layout Standard -116.001290834057 -115.250965599594 8.095072332381373E-005 \end_layout \begin_layout Standard 169.725277983610 8 4 0 \end_layout \begin_layout Standard \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset \end_layout \begin_layout Standard At the end of the simulation, the final configurations are written into external files (in our case ' \shape italic conf_00?? \shape default ' where '??' marks the replica index. Further information for re-starts are written into ' \shape italic par_R.in \shape default ' that also collects statistics on the run. \end_layout \begin_layout Subsection Parallel tempering on a serial machine. \begin_inset LatexCommand label name "sub:Parallel-tempering-on" \end_inset \end_layout \begin_layout Standard Parallel tempering is useful even if you have only a single processor machine available. The algorithm provides better convergence than canonical Monte Carlo and avoids the scheduling problems of simulated annealing. In this example (see \emph on parallel_tempering_s.f \emph default in the \emph on EXAMPLES \emph default sub directory), we perform a parallel tempering simulation on Met-enkaphalin using 5 replicas. Each replica runs at a different temperature. The temperatures are read in from the file \emph on temperatures \emph default . The values are \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard 1 1000 \end_layout \begin_layout Standard 2 600 \end_layout \begin_layout Standard 3 400 \end_layout \begin_layout Standard 4 300 \end_layout \begin_layout Standard 5 250. \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset We start the simulation from randomized configurations with 100 sweeps for equilibration before we exchange any replicas. The main simulation takes 10,000 sweeps with measurements taken every 10 sweeps. An exchange of replicas is also attempted every 10 sweeps. This gives us the following variables \end_layout \begin_layout Standard \begin_inset ERT status open \begin_layout Standard \backslash begin{verbatim} \end_layout \begin_layout Standard num_rep = 5 \end_layout \begin_layout Standard nequi = 100 \end_layout \begin_layout Standard nsweep = 10000 \end_layout \begin_layout Standard nmes = 1 \end_layout \begin_layout Standard newsta = .true. \end_layout \begin_layout Standard switch = 1 \end_layout \begin_layout Standard \backslash end{verbatim} \end_layout \end_inset and we call the simulation with \begin_inset listings lstparams "breaklines=true,float=h,language=Fortran,showstringspaces=false" inline false status collapsed \begin_layout Standard call partem_s(num_rep, nequi, nsweep, nmes, nmes, newsta, switch) \end_layout \end_inset The data is written into file \emph on time.d \emph default . The final configurations are stored in \emph on conf_000?.var \emph default , where the question mark is replaced by the replica ID. \end_layout \begin_layout Standard \begin_inset LatexCommand bibtex options "unsrt" bibfiles "/home/meinke/research/writeups/bibliography" \end_inset \end_layout \begin_layout Standard \begin_inset LatexCommand printindex \end_inset \end_layout \end_body \end_document