#--------------------------------------------------------------------------# # LICENSE INFO: # #--------------------------------------------------------------------------# # This file is part of CAMPARI. # # # # Version 5.0 # # # # Copyright (C) 2024, The CAMPARI development team (current and former # # contributors) # # Andreas Vitalis, Adam Steffen, Rohit Pappu, Hoang # # Tran, Albert Mao, Xiaoling Wang, Jose Pulido, # # Nicholas Lyle, Nicolas Bloechliger, Marco Bacci, # # Davide Garolini, Jiri Vymetal, Cassiano Langini, # # Alessia Del Panta, Francesco Cocina, Julian # # Widmer # # # # Website: http://sourceforge.net/projects/campari/ # # # # CAMPARI is free software: you can redistribute it and/or modify # # it under the terms of the GNU General Public License as published by # # the Free Software Foundation, either version 3 of the License, or # # (at your option) any later version. # # # # CAMPARI 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. # # # # You should have received a copy of the GNU General Public License # # along with CAMPARI. If not, see . # #--------------------------------------------------------------------------# # AUTHORSHIP INFO: # #--------------------------------------------------------------------------# # # # MAIN AUTHOR: Andreas Vitalis # # # #--------------------------------------------------------------------------# class PMS attr_accessor :mpi_mode # type of MPI calc attr_accessor :pmode # parallel mode attr_accessor :nsteps # total steps attr_accessor :mcflag # whether MC moves attr_accessor :cocho # type of container attr_accessor :dofflag # what degrees of freedom attr_accessor :simtype # type of sampler attr_accessor :holomode # level of holonomic constraints attr_accessor :prepmode # level of cluster data preprocessing attr_accessor :pigsflag # whether PIGS run attr_accessor :dockrun # whether MOL2 docking run attr_accessor :docksrun # whether MOL2 docking run in scoring mode attr_accessor :dockrefm # whether there is a reference molecule attr_accessor :dckstyle # what kind of docking to run (coarsely) attr_accessor :dockteth # what kind of tethered docking to perform attr_accessor :dockouts # what kind of information to extract from a docking run (coarsely) attr_accessor :netcsnps # total snapshots (points) for clustering et al. attr_accessor :ghosting # whether FEG is on attr_accessor :absinth # whether ABSINTH is on attr_accessor :stdnbs # whether standard set of NB interactions is present attr_accessor :patching # whether patches wanted attr_accessor :talign # whether trajectory alignment is used attr_accessor :havetmpl # whether a use of template is prescribed attr_accessor :emptyh # whether Hamiltonian is left empty attr_accessor :cdistance# distance metric for clustering attr_accessor :hfourier # related to distance metric for clustering attr_accessor :cdimred # dimensionality reduction for clustering attr_accessor :havehsl # whether HSL is linked attr_accessor :posttask # type of post-processing task attr_accessor :ensemble # type of thermodynamic ensemble to use end class String def is_numeric? begin Float(self) rescue false # not numeric else true # numeric end end def is_integer? begin Float(self) rescue false else if (self[0..0] == "-") if [*self[1..-1].match('[0-9]*')].join == self[1..-1] true else false end else if [*self[0..-1].match('[0-9]*')].join == self[0..-1] true else false end end end end end def curated_bin_answer(inputstr) val = -1 while val < 0 if ((inputstr.downcase)[0..0] == "y") || (inputstr[0..0] == "1") val = 1 elsif ((inputstr.downcase)[0..0] == "n") || (inputstr[0..0] == "0") val = 0 else puts "Please enter an allowed value:" inputstr = gets.chomp end end return val end def curated_int_answer(inputstr,miva,mava) val = miva - 1 while (inputstr.is_integer? == false) puts "Please enter a correctly formatted integer:" inputstr = gets.chomp end val = Integer(inputstr) while ((val < miva) || (val > mava)) puts "Please enter an allowed value:" inputstr = gets.chomp while (inputstr.is_integer? == false) puts "Please enter a correctly formatted integer:" inputstr = gets.chomp end val = Integer(inputstr) end return val end def curated_num_answer(inputstr,miva,mava) notdone = true while (notdone == true) if (inputstr.is_numeric? == true) if ((Float(inputstr) > miva) && (Float(inputstr) < mava)) notdone = false return Float(inputstr) else puts "Please enter a number in the permissible range." end else puts "Please enter a correctly formatted number." end inputstr = gets.chomp end return Float("0.0") end def nest1_parallelism(keyf,pmo) puts "\nCAMPARI can utilize CPU resources in parallel in different ways. While the OpenMP layer is largely abstract (requires low user awareness), the MPI layer is strongly linked to the type of task you want CAMPARI to perform. Which of these options do you require?\n(0): No parallelism whatsoever (single copy of a system/task running on a single thread)\n(1): OpenMP (threads)-accelerated execution (a single copy/task runs using multiple threads in a shared memory-environment, e.g., a single desktop computer)\n(2): MPI (distributed memory) parallelism only (multiple copies of a system or task are executed simultaneously using MPI on a single or on multiple computers, but a single copy uses only a single thread)\n(3): Hybrid OpenMP/MPI parallelism (multiple copies of a system or task are executed simultaneously using MPI on a single or on multiple computers, and each copy uses several threads)\n" keyf.puts("#\n# Parallel settings") pmo.mpi_mode = 0 pmo.pmode = curated_int_answer(gets.chomp,0,3) case pmo.pmode when 0 keyf.puts("FMCSC_NRTHREADS 1 # listed here just so that an inadvertent use of the OpenMP executable behaves as much as possible like the serial one") when 1 keyf.puts("# OpenMP settings\n# The 'right' number of threads depends on a lot of things. For expensive calculations, short benchmark runs are recommended.\n# Good numbers in general are those that line up with processor architecture and allow hyperthreading when supported (e.g., one socket with a modern 8-core CPU could use 8 or 16 threads).\n# On workstations, it may be advisable to not occupy all possible hardware threads in order to keep the usability of the machine appropriate.") keyf.puts("FMCSC_NRTHREADS ???? # the number of OpenMP threads used to accelerate the run\nFMCSC_NRTHREADSONMOL 4 # if you have long chains and enough threads, a larger value might be useful); with fewer threads, the value is adjusted automatically\nFMCSC_THREADS_DLB_FREQ 1000 # the number in steps for how frequently to restart dynamic load balancing (DLB)\nFMCSC_THREADS_DLB_STOP 500 # the number in steps for how long to maximally extend a DLB stretch\n# FMCSC_THREADS_DLB_EXT 1 # in rare cases, it may help to estimate load in an average sense across more than 1 step\nFMCSC_THREADS_VERBOSE 1 # increase from 1 to get more output to the OpenMP task distribution to threads\n# FMCSC_THREADS_TEST 0 # a debugging keyword largely for developer use") when 2 keyf.puts("# For extra safety, make sure that the environment variable OMP_NUM_THREADS is 1 so that an inadvertent use of the OpenMP/MPI executable behaves as much as possible like the pure MPI one.") nest2_mpiism(keyf,pmo) when 3 keyf.puts("# OpenMP settings\n# In a hybrid OpenMP/MPI calculation, a fixed number of resources (cores) are assigned to MPI processes and the threads they start.\n# This is a nontrivial problem, and for CAMPARI we generally recommend to pin threads such that an MPI processes occupies a hardware-local group of threads (e.g., all from the same socket).\n# Because of cache limitations, CAMPARI will tend to benefit from hyperthreading, which adds an additional layer of complexity.\n# The exact distribution/pinning is usually handled by a scheduler like SLURM (see tutorial 7 for an example).\n# Importantly, in hybrid mode, the number of threads per MPI process can only be set through the standard environment variable OMP_NUM_THREADS and not through a keyword.") keyf.puts("# FMCSC_NRTHREADS ???? # cannot use keyword: the standard environment variable OMP_NUM_THREADS determines the actual number of threads\nFMCSC_THREADS_DLB_FREQ 1000 # the number in steps for how frequently to restart dynamic load balancing (DLB)\nFMCSC_THREADS_DLB_STOP 500 # the number in steps for how long to maximally extend a DLB stretch\n# FMCSC_THREADS_DLB_EXT 1 # in rare cases, it may help to estimate load in an average sense across more than 1 step\nFMCSC_THREADS_VERBOSE 1 # increase from 1 to get more output to the OpenMP task distribution to threads\n# FMCSC_THREADS_TEST 0 # a debugging keyword largely for developer use") nest2_mpiism(keyf,pmo) end puts "Keywords related to the parallel execution mode have been added." end def nest1_container(keyf,pmo) puts "\nEvery native CAMPARI run needs to define a container for the system. This information may not be relevant for all analysis features but it is for some. The four most common options are:\n(1): Rectangular box with 3D periodic boundary conditions\n(2): Spherical droplet with a restraining wall acting on all atoms\n(3): A 1D periodic cylinder (along the axis) with restraining walls acting on all atoms elsewhere.\n(4): A 3D-periodic triclinic unit cell.\nPlease select one (you can customize less common options manually from the existing keywords)?" choice = curated_int_answer(gets.chomp,1,4) pmo.cocho = choice keyf.puts("#\n# The choice of container is critical in most situations: most interatomic distances are susceptible to periodicity corrections (minimum image).\n# Also, soft-wall potentials create what is effectively an external potential (depends on absolute coordinates).") case choice when 1 keyf.puts("FMCSC_SHAPE 1 # rectangular box\nFMCSC_BOUNDARY 1 # full 3D periodic boundary conditions\nFMCSC_SIZE ???? ???? ???? # the box lengths in xyz in Angstrom\n# FMCSC_ORIGIN ???? ???? ???? # the absolute position of the corner of the box with minimal values for xyz in Angstrom; in full PBC usually not relevant") when 2 keyf.puts("FMCSC_SHAPE 2 # sphere\nFMCSC_BOUNDARY 4 # a half-harmonic restraining potential will keep every atom in the sphere\nFMCSC_SOFTWALL 0.1 # the strength of this potential in kcal/mol/Angstrom^2 ; this is a relatively small value which means that the true volume will be larger than the formal one\nFMCSC_SIZE ???? # the sphere's radius in Angstrom\nFMCSC_ORIGIN ???? ???? ???? # the absolute position of the center of the sphere in Angstrom (xyz)") when 3 keyf.puts("FMCSC_SHAPE 3 # cylinder\nFMCSC_BOUNDARY 1 # the cylinder is periodic along the z-dimension (its axis) and a half-harmonic restraining potential will keep every atom radially restrained\nFMCSC_SOFTWALL 0.1 # the strength of this potential in kcal/mol/Angstrom^2 ; this is a relatively small value which means that the true volume will be larger than the formal one\nFMCSC_SIZE ???? ???? # the cylinder's radius and height in Angstrom\nFMCSC_ORIGIN ???? ???? ???? # the absolute position of the center of the central circular cross section in Angstrom (xyz)") when 4 keyf.puts("FMCSC_SHAPE 4 # triclinic\nFMCSC_BOUNDARY 1 # 3D periodic boundary conditions\nFMCSC_BOXVECTOR1 ???? ???? ???? # the first box vector in Angstrom, often picked by convention as l_x 0.0 0.0\nFMCSC_BOXVECTOR2 ???? ???? ???? # the second box vector in Angstrom, see documentation for how to construct typical shapes like truncated octahedra or rhombic dodecahedra\nFMCSC_BOXVECTOR3 ???? ???? ???? # the third box vector in Angstrom\nFMCSC_ORIGIN ???? ???? ???? # the absolute position of the center of the 0,0,0 corner in Angstrom (xyz)\n# FMCSC_BOXROTATE 1 # use this to force alignment so that the first box vector is along the x-axis on output") end puts "Keywords related to the simulation container have been added.\n\n" end def nest1_mol2mode(keyf,pmo) keyf.puts("#\n# Settings specific to the small molecule screen (and some general settings required in this context).\n# Make sure your sequence file terminates with XXX_N_C where XXX is any 3-letter residue name not corresponding to a CAMPARI-supported residue\nFMCSC_MOL2SCREEN 1 # mandatory flag to turn on this execution mode\nFMCSC_UNSAFE 1 # pretty much universally required in these types of runs\nFMCSC_MOL2FILE ???? # mandatory input file with a library of (small) molecules in mol2-format\nFMCSC_MOL2ISSPLIT 0 # 1 would mean CAMPARI expects MPI-style, separate input files per replica\nFMCSC_SYBYLLJMAP ???? # it is strongly recommended to use abs4.2.ljmap as shipped with CAMPARI in conjunction with abs4.2* parameters\n#\nFMCSC_MOL2CUTMODE 2 # define residue radii (for cutoff screening) independent of conformation\nFMCSC_MOL2MAXSIZE 500 # maximum allowed number of atoms in a small molecule\nFMCSC_MOL2MAXTIME ???? # set an execution time limit in h to create a clean termination (should be smaller than a true termination limit because it will need some buffer)\n#\nFMCSC_MOL2AUXINDEX automatic # special mode for requesting that moving receptor coordinates be written to mol2\n# FMCSC_MOL2AUXINPUT ???? # possibly necessary index file to deal with inputs containing also receptor coordinates") puts "\nThe small molecule screen for simulation-based docking works by repeating a base calculation for every entry (small molecule) in an input file in mol2-format. Thus, it creates a partially dynamic system representation and deals with the vast majority of setup tasks for the dynamic part automatically.\n" pmo.dockrun = 1 if (((pmo.pmode == 0) || (pmo.pmode == 2)) && (pmo.mpi_mode == 2)) puts "However, as a special variant, it can also be used to 'score' pre-existing small molecule conformations, which is essentially an elaborate single-point energy calculation. Do you want to perform a sampling (1) or a scoring (2) run?\n" choice = curated_int_answer(gets.chomp,1,2) case choice when 1 keyf.puts("FMCSC_MOL2SCOREONLY 0 # deactivates scoring mode (default)\n") when 2 pmo.docksrun = 1 pmo.dockrun = 0 keyf.puts("FMCSC_MOL2SCOREONLY 1 # activates scoring mode, which makes most other functions irrelevant\n# FMCSC_MOL2_PDB_RELAXED ???? # static PDB for referencing (only relevant if dynamic receptor coordinates are present)\nFMCSC_NRSTEPS 1 # no steps needed\nFMCSC_DISABLE_ANALYSIS 1 # no analyses available\nFMCSC_MOL2VERBOSE 2 # reasonable level of reporting\nFMCSC_MCCUTMODE 2 # residue-based truncation (this tends to create a negative offset relative to choice 1)\nFMCSC_MOL2RESPECTBOND 3 # always maximize stability of parameter assignments in scoring runs\nFMCSC_N2LOOP 0 # do not use N^2 energy as MC reference energy (primarily for cost)") end else puts "Selected a sampling run for the small molecule screen (hint: to have access to the scoring mode, select a threads-free\n# execution mode in conjunction with the 'replica exchange'-like option for parallelism (should MPI be requested)." end if (pmo.dockrun == 1) if (pmo.mpi_mode == 2) keyf.puts("# Using FMCSC_REMC in conjunction with FMCSC_MOL2SCREEN leads to a distributed MPI screen, i.e., the master replica is in\n# charge of distributing the molecules found in FMCSC_MOL2FILE to the remaining replicas.") elsif (pmo.mpi_mode == 1) keyf.puts("# Using FMCSC_MPIAVG in conjunction with FMCSC_MOL2SCREEN leads to a small molecule screen where the MPI layer is embedded\n# which means that all replicas process all molecules in coordinated fashion (including data-pooling).") end puts "For how many total elementary steps do you want to run an individual simulation for an individual molecule (simulation settings are set further along)?\n" pmo.nsteps = curated_int_answer(gets.chomp,1,1000000000) keyf.print("FMCSC_NRSTEPS ",pmo.nsteps," # this will be a difficult trade-off to strike for large numbers of screened molecules\nFMCSC_EQUIL ",(pmo.nsteps/10).to_i," # depending on how data are collected this might not be the best choice\nFMCSC_DISABLE_ANALYSIS 1 # since most analyses are not supported, this is the best starting point\n") puts "Do you have a reference molecule (also in mol2-format) available that you want to use to guide the docking?\n" if (curated_bin_answer(gets.chomp) == 1) pmo.dockrefm = 1 keyf.puts("FMCSC_MOL2_REFMOL ???? # the input file with the reference molecule\nFMCSC_MOL2AUXINDEX automatic # special mode for outputting in mol2") else pmo.dockrefm = 0 keyf.puts("FMCSC_MOL2AUXINDEX automatic # special mode") end puts "Which description best matches your application?\n 1) I want to do rigid docking of a small molecule library to a well-defined binding site.\n 2) I want to do flexible docking of a small molecule library to a well-defined binding site.\n 3) I want to run a normal simulation with one or few small molecules.\n 4) I have just small molecules, for example, for conformer generation.\n 5) Something else.\n" pmo.dckstyle = curated_int_answer(gets.chomp,1,5) case pmo.dckstyle when 1 keyf.puts("FMCSC_MOL2FRZMODE 3 # only rigid-body coordinates can change in small molecules\nFMCSC_MOL2RANDOMIZE 1 # randomizes rigid-body coordinates only when no tethering occurs\nFMCSC_MOL2RESPECTBOND 0 # using 1-3 is usually better but it depends on the input file\nFMCSC_MOL2FOCUS -0.1 # freezes everything on the receptor side") when 2 keyf.puts("FMCSC_MOL2FRZMODE 2 # everything is left free in small molecules\nFMCSC_MOL2RANDOMIZE 1 # randomizes rigid-body coordinates only when no tethering occurs, you might consider switching to 4\nFMCSC_MOL2RESPECTBOND 0 # using 1-3 is usually better but it depends on the input file") case pmo.dockrefm when 1 keyf.puts("FMCSC_MOL2FOCUS ???? # freezes everything on the receptor side") when 0 keyf.puts("# FMCSC_MOL2FOCUS ???? # this would be the solution with a well-defined reference molecule\nFMCSC_FRZFILE ???? # put constraints, to leave very specific degrees of freedom free to move, use 'A' mode") end when 3 keyf.puts("FMCSC_MOL2FRZMODE 2 # everything is left free in small molecules\nFMCSC_MOL2RANDOMIZE 0 # keep full input information\nFMCSC_MOL2RESPECTBOND 1 # 1 is a good option for well-curated molecules, use 3 if you start from prior docking poses generated by CAMPARI\nFMCSC_MOL2FOCUS 1.0e10 # leave everything free to move on the receptor side\nFMCSC_FRZFILE ???? # desired constraints should be explicit") when 4 keyf.puts("FMCSC_MOL2FRZMODE 2 # everything is left free in small molecules\nFMCSC_MOL2RANDOMIZE 2 # randomize rigid-body\nFMCSC_MOL2DRESTMODE 0 # nothing to restrain\nFMCSC_MOL2RESPECTBOND 0 # if this is ligand sampling to create finite differences, definitely use 3") when 5 keyf.puts("FMCSC_MOL2FRZMODE 0 # defer to underlying settings\nFMCSC_MOL2RANDOMIZE 0 # keep full input information\nFMCSC_MOL2RESPECTBOND 1 # if this is ligand sampling to create finite differences, definitely use 3\n") end if (pmo.dckstyle != 4) puts "If the reference molecule or the first molecule in the library file is tagged to indicate a reference substructure, how do you want to use this information?\n 1) There is no tagging.\n 2) I want only exact matches for substructure alignment and tethered docking.\n 3) I want to use alignment primarily as a way to place molecules in the binding site.\n 4) I want unusual matches to be explored, and I realize that this can slow down my calculation drastically.\n 5) Something else." pmo.dockteth = curated_int_answer(gets.chomp,1,4) case pmo.dockteth when 1 if (pmo.dckstyle <= 2) keyf.puts("# FMCSC_SC_DREST 0.1 # consider setting up some mild restraints\n# FMCSC_MOL2DRESTMODE 2 # without tagging, they have to be taken from heavy atom positions (envelope)\nFMCSC_MOL2DRESTBUF 5.0 # buffer in Angstrom") end when 2 keyf.puts("FMCSC_SC_DREST 0.1 # global scale factor for position restraints\nFMCSC_MOL2DRESTMODE 1 # use the tagged atoms to create specific restraints\nFMCSC_MOL2DRESTBUF 1.0 # allow a bit of leeway (value in Angstrom)\nFMCSC_MOL2ASSIMILAR 0.25 # be stringent on chemical matches\nFMCSC_MOL2RMSDTHRESH 0.5 # be stringent on geometric matches") when 3 keyf.puts("FMCSC_SC_DREST 0.1 # global scale factor for position restraints\nFMCSC_MOL2DRESTMODE 2 # use all heavy atoms to create envelope-type restraints\nFMCSC_MOL2DRESTBUF 5.0 # allow leeway (value in Angstrom)\nFMCSC_MOL2ASSIMILAR 2.5 # be lenient on chemical matches\nFMCSC_MOL2RMSDTHRESH 0.5 # be stringent on geometric matches to avoid duplicities") when 4 keyf.puts("FMCSC_SC_DREST 0.1 # global scale factor for position restraints\nFMCSC_MOL2DRESTMODE 1 # use the tagged atoms to create specific restraints\nFMCSC_MOL2DRESTBUF 2.0 # allow leeway (value in Angstrom)\nFMCSC_MOL2ASSIMILAR 3.5 # be lenient on chemical matches\nFMCSC_MOL2RMSDTHRESH 2.5 # be lenient on geometric matches (do not include symmetric hydrogens in tagged substructure") when 5 if (pmo.dckstyle <= 2) keyf.puts("# FMCSC_SC_DREST ???? # consider setting up some mild restraints\n# FMCSC_MOL2DRESTMODE ???? # without tagging, they have to be taken from heavy atom positions (envelope), with tagging they can come from tagged atoms\nFMCSC_MOL2DRESTBUF ???? # buffer in Angstrom") end end end puts "And which description best matches your desired outputs?\n 1) I just want best poses.\n 2) I want full disclosure, and I realize that this might create an enormous amount of data very quickly.\n 3) I want to save docking-specific information with a focus on pose diversity.\n 4) Something else.\n" pmo.dockouts = curated_int_answer(gets.chomp,1,4) case pmo.dockouts when 1 keyf.puts("FMCSC_MOL2OUTMODE 3 # only minimum energy poses stored\nFMCSC_CCOLLECT 2000000000 # do not collect data for pose clustering\nFMCSC_MOL2VERBOSE 1 # only basic reporting\nFMCSC_MOL2PDBMODE 0 # do not save *END.pdbs per molecule") when 2 keyf.print("FMCSC_MOL2OUTMODE 2 # final, minimum energy, and cluster centroid poses are stored\nFMCSC_MOL2ENMODE 2 # averaged energies are reported as scores: use a scoring run to get single-point values for all stored poses\nFMCSC_MOL2CLUMODE 2 # a reduced set will usually suffice, irrespective of the distance metric\nFMCSC_CCOLLECT ",(pmo.nsteps/1000).to_i," # this a guess meant to produce 900 analyzed instantaneous poses\nFMCSC_CFILE ???? # you will have to create a dummy file\nFMCSC_CDISTANCE 5 # RMSD might be the simplest but only if your receptor remains in place, otherwise change to 2 ot 7\nFMCSC_CALIGN 0 # receptor must remain fixed\nFMCSC_CRADIUS 2.0 # too fine a clustering will be a bad idea\nFMCSC_MOL2THRESH 0.1 # 10% threshold usually means up to 3 clusters are retained\nFMCSC_MOL2VERBOSE 4 # maximum level of reporting\nFMCSC_MOL2PDBMODE 3 # save *END.pdbs and, possibly, *RELAXED.pdbs per molecule\n") when 3 keyf.print("FMCSC_MOL2OUTMODE 1 # minimum energy and cluster centroid poses are stored\nFMCSC_MOL2ENMODE 2 # averaged energies are reported as scores\nFMCSC_MOL2CLUMODE 2 # a reduced set will usually suffice, irrespective of the distance metric\nFMCSC_CCOLLECT ",(pmo.nsteps/1000).to_i," # this a guess meant to produce 900 analyzed instantaneous poses\nFMCSC_CFILE ???? # you will have to create a dummy file\nFMCSC_CDISTANCE 5 # RMSD might be the simplest but only if your receptor remains in place, otherwise change to 2 ot 7\nFMCSC_CALIGN 0 # receptor must remain fixed\nFMCSC_CRADIUS 2.0 # too fine a clustering will be a bad idea\nFMCSC_MOL2THRESH 0.05 # 5% threshold usually means up to 5 clusters are retained\nFMCSC_MOL2VERBOSE 2 # reasonable level of reporting\nFMCSC_MOL2PDBMODE 0 # do not save *END.pdbs per molecule\n") when 4 keyf.print("FMCSC_MOL2OUTMODE ???? # 1 is the most common option, 4 is rarely used (trajectory-like output)\nFMCSC_MOL2ENMODE ???? # can either use single-point or averaged energies as scores\n# FMCSC_MOL2CLUMODE 2 # a reduced set will usually suffice, irrespective of the distance metric\n# FMCSC_CCOLLECT ",(pmo.nsteps/1000).to_i," # this a guess meant to produce 900 analyzed instantaneous poses\n# FMCSC_CFILE ???? # you would have to create a dummy file\n# FMCSC_CDISTANCE 5 # RMSD might be the simplest but only if your receptor remains in place, otherwise change to 2 ot 7\n# FMCSC_CALIGN 0 # receptor must remain fixed\n# FMCSC_CRADIUS 2.0 # too fine a clustering will be a bad idea\n# FMCSC_MOL2THRESH 0.05 # 5% threshold would usually mean that up to 5 clusters are retained\nFMCSC_MOL2VERBOSE ???? # 0-4 is available\nFMCSC_MOL2PDBMODE 0 # do not save any pdb per molecule\n") end end puts "Keywords related to the small molecule screen have been added.\n\n" end def nest1_stype(keyf,pmo) if (pmo.mpi_mode == 0) puts "\nSince this is not an analysis run, what type of simulation do you want to perform?" else puts "\nSince this is not an analysis run, what type of simulation do you want an individual replica (MPI process) to perform?" end keyf.puts("#\n# The settings below define the type of simulation (sampler, bath coupling, etc).") pmo.mcflag = 0 pmo.holomode = 0 terms_unselected = Array.new(7){ |i| i+1 } terms_unselected.delete(4) if (pmo.mpi_mode > 0) terms_unselected.delete(6) end choice = 0 terms_unselected.each do |i| printf("(%2d): ",i) case i when 1 puts "Pure Monte Carlo (MC) sampling" when 2 puts "Pure Newtonian dynamics sampling" when 3 puts "Pure stochastic dynamics (Langevin) sampling" when 4 puts "Pure Brownian dynamics sampling" when 5 puts "Hybrid Newtonian dynamics/MC sampling" when 6 puts "Minimization" when 7 puts "Hybrid stochastic dynamics/MC sampling" end end while (choice == 0) choice = curated_int_answer(gets.chomp,1,7) while (!terms_unselected.include?(choice)) puts "Please select an available option." choice = curated_int_answer(gets.chomp,1,7) end end pmo.simtype = choice keyf.puts("# Simulation extent settings\n") if (pmo.dockrun != 1) puts "Is this a simulation you will need to restart later on?\n" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_RESTART 0 # setting this to 1 will later allow you to restart (extend) a simulation\nFMCSC_RST_MC2MD 0 # setting this to 1 can help with restarting simulations with other samplers") end puts "How many total elementary steps should the simulation contain? If this is a simulation you will need to restart continuously, just enter the value for the first increment.\n" pmo.nsteps = curated_int_answer(gets.chomp,1,2000000000) keyf.print("FMCSC_RSTOUT ",(pmo.nsteps/100).to_i," # for very fast simulations, this should be increased further\nFMCSC_NRSTEPS ",pmo.nsteps," # the fundamental target simulation extent in elementary steps\nFMCSC_EQUIL ???? # the number of elementary steps to truncate at the beginning for analysis tasks and trajectory writing\n# The type of sampler\n") else keyf.print("FMCSC_RSTOUT 1000000000 # restarts are not possible in small molecule screens\n# The type of sampler\n") end case choice when 1 keyf.puts("FMCSC_DYNAMICS 1 # 1 selects pure Monte Carlo: this implies internal coordinates\nFMCSC_MCCUTMODE 1 # consider changing this to 2 to switch to residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # you can enable this in pure MC runs but the cutoff-free energies might still be misleading\nFMCSC_USESCREEN 0 # using the screen requires good knowledge of how it works, 0 disables it\n# FMCSC_BARRIER ???? # required if you use the screen") pmo.mcflag = 1 nest2_moveset(keyf,pmo) pmo.dofflag = 1 if (pmo.mpi_mode != 2) nest3_wl(keyf,pmo) end when 2 if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) keyf.puts("FMCSC_DYNAMICS 2 # 2 selects ballistic (Newtonian) dynamics\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # never compute N^2 energy") else keyf.puts("FMCSC_DYNAMICS 2 # 2 selects ballistic (Newtonian) dynamics\nFMCSC_SYSFRZ 2 # remove translational drift from the system (can avoid equipartition problems)\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # never compute N^2 energy") end nest2_ensemble(keyf,pmo) nest2_dofs(keyf,pmo,1) nest2_relax(keyf,pmo) nest2_coupling(keyf,pmo) when 3 keyf.puts("FMCSC_DYNAMICS 3 # 3 selects stochastic dynamics in low friction regime\nFMCSC_SYSFRZ 1 # the removal of translational drift from the system is usually unnecessary in stochastic dynamics\nFMCSC_ENSEMBLE 1 # stochastic dynamics imply temperature bath coupling\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # never compute N^2 energy") nest2_dofs(keyf,pmo,1) nest2_relax(keyf,pmo) nest2_coupling(keyf,pmo) when 5 pmo.mcflag = 1 if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) keyf.puts("FMCSC_DYNAMICS 5 # 5 selects mixed Monte Carlo / Newtonian dynamics sampling\nFMCSC_ENSEMBLE 1 # the MC part implies temperature bath coupling\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # do not use N^2 energy as MC reference energy") else keyf.puts("FMCSC_DYNAMICS 5 # 5 selects mixed Monte Carlo / Newtonian dynamics sampling\nFMCSC_SYSFRZ 2 # the removal of translational drift is technically incompatible with MC but can avoid equipartition problems\nFMCSC_ENSEMBLE 1 # the MC part implies temperature bath coupling\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # do not use N^2 energy as MC reference energy") end pmo.dofflag = 1 nest2_relax(keyf,pmo) nest2_coupling(keyf,pmo) nest2_hybrids(keyf,pmo) nest2_moveset(keyf,pmo) when 6 keyf.puts("FMCSC_DYNAMICS 6 # 6 selects minimization\n# FMCSC_ENSEMBLE 1 # if at all in an ensemble, CAMPARI minimizers operate in the NVT one\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # never compute N^2 energy") nest2_dofs(keyf,pmo,0) nest2_relax(keyf,pmo) nest2_mini(keyf,pmo) when 7 pmo.mcflag = 1 keyf.puts("FMCSC_DYNAMICS 7 # 7 selects mixed Monte Carlo / stochastic dynamics in low friction regime\nFMCSC_SYSFRZ 1 # the removal of translational drift from the system is usually unnecessary in stochastic dynamics\nFMCSC_ENSEMBLE 1 # stochastic dynamics and MC both imply temperature bath coupling\nFMCSC_MCCUTMODE 2 # residue-based truncation (the generally supported case)\nFMCSC_N2LOOP 0 # do not use N^2 energy as MC reference energy\n") pmo.dofflag = 1 nest2_relax(keyf,pmo) nest2_hybrids(keyf,pmo) nest2_moveset(keyf,pmo) nest2_coupling(keyf,pmo) end if ((pmo.dockrun == 1) && (pmo.dckstyle != 3) && (pmo.simtype != 1)) keyf.puts("FMCSC_THRESHOLD_INCR 0.1 # make sure unstable simulations can terminate without causing unhandled exceptions\nFMCSC_THRESHOLD_TEMP 5000.0 # another protection mechanism") elsif (pmo.simtype != 1) keyf.puts("# FMCSC_THRESHOLD_TEMP 600.0 # if you are solvating a macromolecule with a dense fluid or are starting from a clashy structure in general, this can be very useful in dynamics-based methods to keep the run going: definitely disable it once it is no longer needed\n") end if ((pmo.dockrun == 1) && ((pmo.dckstyle <= 2) || (pmo.dckstyle == 4))) keyf.puts("FMCSC_SKIPFRZ 1 # if it is supported, great, if not it will be disabled anyway") puts "Keywords related to the choice of simulation/sampler have been added.\n\n" else nest3_frzs(keyf,pmo) puts "Keywords related to the choice of simulation/sampler have been added.\n\n" end end def nest1_afiletype(keyf,pmo) keyf.puts("FMCSC_PDB_READMODE 2 # uses coordinate information directly and assumes all atoms present (implied default in trajectory analysis mode)") puts "Do you want to analyze (1): pdb trajectory (MODEL/ENDMDL); (2): pdb-files (systematically numbered); (3): xtc trajectory; (4): dcd trajectory; (5): NetCDF trajectory; (6): PSQL database entry?" choice = curated_int_answer(gets.chomp,1,6) case choice when 1 keyf.puts("FMCSC_PDB_FORMAT 1 # 1 selects a PDB trajectory\nFMCSC_PDBFILE ???? # the pdb trajectory file\nFMCSC_PDB_R_CONV ???? # you need to change from the default of 1 only for specific issues: 3 is good when you work with files directly from CHARMM, 4 with files directly from AMBER") when 2 keyf.puts("FMCSC_PDB_FORMAT 2 # 2 selects individual PDB files\nFMCSC_PDBFILE ???? # should be first of systematically numbered pdb files\nFMCSC_PDB_R_CONV ???? # you need to change from the default of 1 only for specific issues: 3 is good when you work with files directly from CHARMM, 4 with files directly from AMBER") when 3 keyf.puts("FMCSC_PDB_FORMAT 3 # 3 selects a CAMPARI/GROMACS format xtc-trajectory\nFMCSC_XTCFILE ???? # the .xtc trajectory file; binary files need a template unless the order of atoms is exactly the same as CAMPARI default") when 4 keyf.puts("FMCSC_PDB_FORMAT 4 # 4 selects a CAMPARI/CHARMM/NAMD format dcd-trajectory\nFMCSC_DCDFILE ???? # the .dcd trajectory file; binary files need a template unless the order of atoms is exactly the same as CAMPARI default") when 5 keyf.puts("FMCSC_PDB_FORMAT 5 # 5 selects an CAMPARI/AMBER-format NetCDF trajectory\nFMCSC_NETCDFFILE ???? # the NetCDF (.nc) trajectory file; binary files need a template unless the order of atoms is exactly the same as CAMPARI default") when 6 keyf.puts("FMCSC_PDB_FORMAT 6 # 6 selects the PSQL database format and entry\nFMCSC_PSQL_R_TRAJKEY ???? # name of the targeted data set (in the so-called ensemble tier) in the database\nFMCSC_PSQL_R_DATABASE ???? # the name of the database (full read permissions required)\nFMCSC_PSQL_R_DBHOST ???? # the database host's (IP) address, default is 127.0.0.1 (local)\n# FMCSC_PSQL_R_DBPORT 5432 # change only if your database administrators changed the default port\n# FMCSC_PSQL_R_SIMUOFF ???? # if you don't want to analyze all trajectories associated with the key but a specific one\nFMCSC_PSQL_R_RANK 0 # this is the default and the only choice that allows parallel analyses on different trajectories: change only with an exact purpose in mind\n# FMCSC_PSQL_R_LIMIT ???? # you could limit the number of trajectories obtained from the database, rarely needed\n# FMCSC_PSQL_R_PREFSORT 2 # this is the default order (rank first) and the more intuitive one") end keyf.puts("FMCSC_NRSTEPS ???? # the number of the snapshots in the input trajectory (or the number of pdb files) should be exact (if it is unknown choose a safe larger value)\nFMCSC_EQUIL 0 # change from 0 if you want to discard snapshots at the beginning\n# FMCSC_FRAMESFILE ???? # uncomment and provide a valid file if you want to select specific snapshots for analysis or apply sampling weights") puts "Keywords related to the choice of input file have been added.\n\n" if (choice > 2) keyf.puts("# If you have a binary format and unsupported residues, FMCSC_PDB_TEMPLATE might have to serve two purposes from the same file.\n# In these cases, it is often cleaner to use FMCSC_PDB_ATOMMAP for the remapping function.") keyf.puts("FMCSC_PDB_TEMPLATE ???? # the template is one of two ways to remap the order of atoms from binary trajectory files written by other software\n# FMCSC_PDB_ATOMMAP ???? # this input file can replace the function of the template for remapping and is necessary if, for example, atoms in the trajectory file are not grouped by residue as CAMPARI expects it\nFMCSC_PDB_R_CONV ???? # you need to change from the default of 1 only for specific issues: 3 is good when you work with files directly from CHARMM, 4 with files directly from AMBER\n# FMCSC_PDB_INPUTSTRING \"a6,a5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3\" # you need to activate this if your input file is a non-standard PDB") pmo.havetmpl = 1 else puts "Does the system contain residues not supported by CAMPARI (y/n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_PDB_TEMPLATE ???? # a template is always needed to infer the molecular topology of unsupported residues") pmo.havetmpl = 1 end end if (pmo.cocho != 2) puts "Are the molecules in the input trajectories 'broken' such that all atoms always reside in the central unit cell (GROMACS-style) (y) or are they 'whole' but stick out of the unit cell (n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_XYZ_FORCEBOX 1 # 1 assumes molecules are broken into pieces from multiple images on input") else keyf.puts("FMCSC_XYZ_FORCEBOX 0 # 0 assumes molecules are intact on input (not broken into pieces from multiple images)") end end puts "Is this trajectory obtained at custom volume (n) or do volume changes need to be read from the trajectory data (y)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_ENSEMBLE 3 # assume NPT conditions, currently only in analysis mode") else keyf.puts("FMCSC_ENSEMBLE 1 # assume NVT conditions") end end def nest1_hamiltonian(keyf,pmo,mode) # define some shortcuts ljs = "FMCSC_SC_IPP 1.0 # scale factor for inverse power potential (excluded volume)\nFMCSC_SC_ATTLJ 1.0 # scale factor for dispersion potential (LJ)\nFMCSC_IPPEXP 12 # 12 is also the default (LJ)\nFMCSC_SIGRULE ???? # please look up correct value for your parameter file in documentation\nFMCSC_EPSRULE ???? # please look up correct value for your parameter file in documentation\nFMCSC_HSSCALE 1.0 # can be used to scale sigma parameters globally" absmodel = "#\n# The short-range interaction model controls what interactions are excluded because of molecular topology\nFMCSC_INTERMODEL 1 # strong nonbonded exclusion model (as used with ABSINTH)\nFMCSC_ELECMODEL 2 # strong short-range exclusion model for electrostatics\nFMCSC_FUDGE_ST_14 1.0 # 14-classified interactions are counted fully for LJ/WCA\nFMCSC_FUDGE_EL_14 1.0 # 14-classified interactions are counted fully for Coulomb\nFMCSC_MODE_14 1 # choice matters only if FMCSC_FUDGE_ST_14 or FMCSC_FUDGE_EL_14 deviate from 1.0" absmodel2 = "#\n# The short-range interaction model controls what interactions are excluded because of molecular topology\nFMCSC_INTERMODEL 1 # strong nonbonded exclusion model (as used with ABSINTH)\nFMCSC_ELECMODEL 3 # strong short-range exclusion model for electrostatics with self-correction (largely the same as 2)\nFMCSC_FUDGE_ST_14 1.0 # 14-classified interactions are counted fully for LJ/WCA\nFMCSC_FUDGE_EL_14 1.0 # 14-classified interactions are counted fully for Coulomb\nFMCSC_MODE_14 1 # choice matters only if FMCSC_FUDGE_ST_14 or FMCSC_FUDGE_EL_14 deviate from 1.0" ffmodel = "#\n# The short-range interaction model controls what interactions are excluded because of molecular topology\nFMCSC_INTERMODEL 2 # weak nonbonded exclusion model (as used with standard force fields)\nFMCSC_ELECMODEL 1 # weak short-range exclusion model for electrostatics\nFMCSC_FUDGE_ST_14 ???? # please look up the correct scale factor for 14-classified interactions (LJ/WCA) for your parameter file in the docu.\nFMCSC_FUDGE_EL_14 ???? # the same for Coulomb interactions\nFMCSC_MODE_14 1 # choice matters only if FMCSC_FUDGE_ST_14 or FMCSC_FUDGE_EL_14 deviate from 1.0" ffmodel2 = "#\n# The short-range interaction model controls what interactions are excluded because of molecular topology\nFMCSC_INTERMODEL 3 # weak nonbonded exclusion model (as used with GROMOS force fields)\nFMCSC_ELECMODEL 1 # weak short-range exclusion model for electrostatics\nFMCSC_FUDGE_ST_14 ???? # please look up the correct scale factor for 14-classified interactions (LJ/WCA) for your parameter file in the docu.\nFMCSC_FUDGE_EL_14 ???? # the same for Coulomb interactions\nFMCSC_MODE_14 1 # choice matters only if FMCSC_FUDGE_ST_14 or FMCSC_FUDGE_EL_14 deviate from 1.0" modeldebug = "# FMCSC_INTERREPORT 0 # you can enable this report for small systems to understand the general short-range interaction model\n# FMCSC_ELECREPORT 0 # you can enable this report for small systems to understand the Coulomb short-range interaction model" pmo.emptyh = 0 puts "\nNext, I would like to know whether any energies/forces should be evaluated at all (y/n)?\nChoosing 'n' will create a minimal set of keywords with no active energy terms.\n" if (curated_bin_answer(gets.chomp) == 0) keyf.puts("#\n# Empty Hamiltonian\nFMCSC_SC_IPP 0.0 # disable only energy term enabled by default\nFMCSC_CUTOFFMODE 4 # there are technical reasons for choosing this") puts "The relevant keywords have been added.\n\n" pmo.ghosting = 0 pmo.absinth = 0 pmo.emptyh = 1 return end puts "Before starting with energy terms, I need to know what short-range interaction model you intend to use for nonbonded interactions. If you do not intend to use any nonbonded, pairwise interactions whatsoever, choose '5'.\n(1): Use the ABSINTH-like short-range interaction model\n(2): Use the ABSINTH-like short-range interaction model with self-correction\n(3): Use the standard force field model\n(4): Use the standard force field model with GROMOS corrections\n(5): There is no need for a nonbonded short-range interaction model\n" choicesr = curated_int_answer(gets.chomp,1,5) case choicesr when 1 keyf.puts(absmodel) keyf.puts(modeldebug) when 2 keyf.puts(absmodel2) keyf.puts(modeldebug) when 3 keyf.puts(ffmodel) keyf.puts(modeldebug) when 4 keyf.puts(ffmodel2) keyf.puts(modeldebug) when 5 # do nothing end keyf.puts("#\n# The truncation model (cutoffs) is the primary determinant of computational speed in the presence of general nonbonded interactions.") puts "Similarly, I need to know what long-range truncation (cutoff) scheme you intend to use for nonbonded interactions. If you do not intend to use any nonbonded, pairwise interactions whatsoever, choose '1'.\n(1): Use topology-based cutoffs (fastest for up to a range of 100s to 1000s residues, normal choice for implicit solvent)\n(2): Use grid-based cutoffs (faster for larger numbers of residues, normal choice for explicit solvent)\n(3): No cutoffs should be applied (do not choose this option unless you know exactly why)\n" choice2 = curated_int_answer(gets.chomp,1,3) case choice2 when 1 keyf.puts("FMCSC_CUTOFFMODE 4 # enables topology-assisted cutoffs of nonbonded interactions\nFMCSC_NBCUTOFF ???? # the general cutoff for LJ/IPP/WCA/SAV; common values for atomic systems are 9-12 Angstrom\nFMCSC_ELCUTOFF ???? # the general cutoff for POLAR/TABUL; common values for atomic systems are 9-14 Angstrom\n") when 2 keyf.puts("FMCSC_CUTOFFMODE 3 # enables grid-based cutoffs of nonbonded interactions\nFMCSC_NBCUTOFF ???? # the general cutoff for LJ/IPP/WCA/SAV; common values for atomic systems are 9-12 Angstrom\nFMCSC_ELCUTOFF ???? # the general cutoff for POLAR/TABUL; common values for atomic systems are 9-14 Angstrom\nFMCSC_GRIDDIM ???? ???? ???? # match to system dimensions so that the cells are cubic\n# FMCSC_GRIDMAXGPNB ???? # you may have to increase this value if you choose a fine grid or long cutoffs (try increasing powers of 2)\n# FMCSC_GRIDMAXRSNB ???? # you may have to increase this value if you choose a very coarse grid\nFMCSC_GRIDREPORT 0 # can be enabled for debugging\n") when 3 keyf.puts("FMCSC_CUTOFFMODE 1 # disables truncation of nonbonded interactions, you should know exactly why you picked this") end case choice2 when 1,2 if ((pmo.simtype != 1) && (pmo.simtype != 6)) keyf.puts("FMCSC_NBL_UP 5 # neighbor list update frequency in dynamics, must be set heuristically (only 1 is rigorously safe but also slowest)") end end puts "I would also like to know whether patching will be required (y/n)? The by far most common scenario requiring patches are unsupported residues (i.e., residues non in the CAMPARI-internal residue database). Patching means to use custom input files to assign or reassign force field parameters." if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) puts "Note that the small molecules are treated automatically and cnanot be patched, so do not answer 'y' unless there are additional unsupported residues (static) that require patches." end pmo.patching = curated_bin_answer(gets.chomp) if (pmo.patching == 1) keyf.puts("# FMCSC_MPATCHFILE ???? # often, mass patches are not really needed") end terms_unselected = Array.new(16){ |i| i+1 } terms_selected = [] puts "A final preparatory question: Do you want to enable the ghosting of (primarily nonbonded) interactions as used in free energy (growth) calculations (y/n)? Choose 'n' if uncertain.\n" pmo.ghosting = curated_bin_answer(gets.chomp) if (pmo.patching == 1) keyf.puts("#\n# Hamiltonian choices\nFMCSC_BIOTYPEPATCHFILE ???? # the best way to patch in CAMPARI is at the biotype level, possibly by means of an appended parameter file\nFMCSC_SC_IPP 0.0 # this is the only energy term that would be enabled by default; setting it to 0.0 explicitly disables it") else keyf.puts("#\n# Hamiltonian choices\n# FMCSC_BIOTYPEPATCHFILE ???? # the best way to patch ff parameters in CAMPARI would be at the biotype level, possibly by means of an appended parameter file\nFMCSC_SC_IPP 0.0 # this is the only energy term that would be enabled by default; setting it to 0.0 explicitly disables it") end if (pmo.ghosting == 1) keyf.puts("FMCSC_GHOST 1 # ghosting of interactions enabled: see special scale factors and custom input file below\nFMCSC_FEGFILE ???? # the required input file for which residues to ghost in the Hamiltonian\nFMCSC_FEG_MODE 2 # relevant only for scaled nonbonded interactions: this is the slightly more versatile option\nFMCSC_FEGREPORT 1 # reports are useful when using custom input files") terms_unselected.delete(3) terms_unselected.delete(5) else keyf.puts("FMCSC_GHOST 0 # no ghosting of (primarily nonbonded) interactions will occur (default)") end if (pmo.dockrun == 1) terms_unselected.delete(14) terms_unselected.delete(11) end choice = 0 pmo.absinth = 0 pmo.stdnbs = 0 needlrel = 0 while (choice != 17) puts "What terms of the energy function should be enabled?" terms_unselected.each do |i| printf("(%2d): ",i) case i when 1 puts "Standard force field terms (includes ABSINTH and vacuum (fully explicit) solvation models)" when 2 puts "Just volume exclusion/dispersion via Lennard-Jones" when 3 puts "Just volume exclusion/dispersion via Weeks-Chandler-Andersen" when 4 puts "Just excluded volume" when 5 puts "Weeks-Chandler-Andersen with Coulomb" when 6 puts "Just excluded volume and Coulomb" when 7 puts "Standard bonded terms (excluding CMAP)" when 8 puts "CMAP terms" when 9 puts "Torsional bias terms (TOR) for individual dihedral angles" when 10 puts "Secondary structure bias terms (ZSEC) based on dihedral angles" when 11 puts "Secondary structure bias terms (DSSP) based on hydrogen bonds" when 12 puts "Polymeric bias terms (POLY) on radii of gyration and asphericities" when 13 puts "Interatomic, tabulated potentials (TABUL) from custom input files" when 14 puts "Distance restraints (DREST) on atom pairs (e.g., NOE-derived restraints)" when 15 puts "Compartmentalization restraints (OSMO): these potentials allow molecules to be confined to subvolumes inside the simulation container" when 16 puts "Spatial density restraints (EMICRO) (e.g., electron microscopy density maps)" end end puts "(17): I do not want to include any of these energy terms or I have selected all the ones I want." choice = curated_int_answer(gets.chomp,1,17) while ((!terms_unselected.include?(choice)) && (choice != 17)) puts "Please select an available option." choice = curated_int_answer(gets.chomp,1,17) end if (choice < 17) terms_selected.push(choice) end if (choice <= 6) (0..6).each do |j| terms_unselected.delete(j) end elsif (choice < 17) terms_unselected.delete(choice) end case choice when 1 case choicesr when 1,2 if (pmo.ghosting == 0) puts "Do you want to use the ABSINTH implicit solvation model (y) or not, i.e., a vacuum background (n)?" if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) puts "(note that the parameterization of ABSINTH for the small molecules is NOT automatic and has to be encoded in the library file)" end if (curated_bin_answer(gets.chomp) == 1) pmo.absinth = 1 end end if (pmo.absinth == 1) keyf.puts("FMCSC_SC_IPP 1.0 # scale factor for inverse power potential (excluded volume)\nFMCSC_SC_ATTLJ 1.0 # scale factor for dispersion potential (LJ)\nFMCSC_IPPEXP 12 # 12 is also the default (LJ)\nFMCSC_SIGRULE 1 # ABSINTH LJ parameters imply this combination rule\nFMCSC_EPSRULE 2 # ABSINTH LJ parameters imply this combination rule\n# FMCSC_HSSCALE 1.0 # can be used to scale sigma parameters globally") if (pmo.patching == 1) keyf.puts("# FMCSC_LJPATCHFILE ???? # in some cases, it may prove useful to selectively patch LJ parameters (rather than biotypes)") end keyf.puts("FMCSC_SC_POLAR 1.0 # scale factor for Coulomb potential\nFMCSC_DIPREPORT 0 # can be used to produce VMD files for visualizing charge group information") if (pmo.patching == 1) keyf.puts("FMCSC_CPATCHFILE ???? # partial charges are often context-dependent, e.g., in molecules with a mix of supported and unsupported residues") else keyf.puts("# FMCSC_CPATCHFILE ???? # partial charges often describe nonlocal (context-dependent) phenomena, e.g., terminal residues in AMBER") end keyf.puts("# The ABSINTH model has many customizable parameters (below) but all deviations from published form(s) must be considered experimental\nFMCSC_SC_IMPSOLV 1.0 # scale factor for the ABSINTH DMFI: any nonzero value also enables charge screening\nFMCSC_SCRMODEL 2 # atom-based screening as published (the group-based model, choice 1, is an alternative)\nFMCSC_IMPDIEL 78.2 # the dielectric can be changed to reflect changes in, for example, temperature\n# FMCSC_FOSMODE 1 # 1 is the default temperature-independent model: can be changed if parameter file is extended to support it\n# FMCSC_FOSREFT ???? # relevant only for temperature-dependent DMFI parameters\nFMCSC_FOSREPORT 0 # you can use it to create a parameter-level report of DMFI parameters\nFMCSC_SAVPROBE 2.5 # the radius of a probe sphere in Angstrom: the solvation shell thickness is twice this value\nFMCSC_SAVMODE 2 # a conservative evolution of the linear approximation of the original model (which would be choice 1)\nFMCSC_FOSFUNC 3 # the original smooth interpolation function: choice 3 is a conservative update with continuous gradients\n# FMCSC_FOSMID ???? # the DMFI mid-point parameter (the published value is the default), largely irrelevant if FMCSC_FOSTAU is large\n# FMCSC_FOSTAU ???? # the DMFI steepness parameter (the published value is the default)\n# FMCSC_FOSTAPER ???? # interval fraction to taper the function over if FOSFUNC is 3, default 0.05\n# FMCSC_FOSSHIFT ???? # related to nonstandard interpolation functions\n# FMCSC_FOSTIGHT ???? # related to nonstandard interpolation functions\n# FMCSC_FOSGRANULE ???? # related to nonstandard interpolation functions\nFMCSC_SCRFUNC 3 # the original smooth interpolation function for charge screening: choice 3 is a conservative update with continuous gradients\n# FMCSC_SCRMID ???? # the charge screening mid-point parameter (the published value is the default), largely irrelevant if FMCSC_SCRTAU is large\n# FMCSC_SCRTAU ???? # the charge screening steepness parameter (the published value is the default)\n# FMCSC_SCRTAPER ???? # interval fraction to taper the function over if SCRFUNC is 3, default 0.05\n# FMCSC_SCRSHIFT ???? # related to nonstandard interpolation functions\n# FMCSC_SCRTIGHT ???? # related to nonstandard interpolation functions\n# FMCSC_SCRGRANULE ???? # related to nonstandard interpolation functions\n# FMCSC_ISQM ???? # related to nonstandard screening models\n# FMCSC_CONTACTDIEL ???? # related to nonstandard screening models\n# FMCSC_SCRMIX ???? # related to nonstandard screening models\n# Note that patching unsupported residues for the ABSINTH model requires a solid grasp of the paradigm.\n# FMCSC_FOSPATCHFILE ???? # this is the file required to set up new or rewritten solvation groups for ABSINTH\n# FMCSC_RPATCHFILE ???? # it is possible to patch atomic radii for SAV/ABSINTH purposes (makes them LJ-independent)\n# FMCSC_SAVPATCHFILE ???? # it is possible to patch maximum SAV fractions per atom if needed\n# FMCSC_ASRPATCHFILE ???? # it is possible to patch atomic reduction factors for SAV calculations") need_lrel = 1 else pmo.stdnbs = 1 keyf.puts(ljs) if (pmo.patching == 1) keyf.puts("# FMCSC_LJPATCHFILE ???? # in some cases, it may prove useful to selectively patch LJ parameters (rather than biotypes)") end keyf.puts("FMCSC_SC_POLAR 1.0 # scale factor for Coulomb potential\nFMCSC_DIPREPORT 0 # can be used to produce VMD files for visualizing charge group information\n# FMCSC_AMIDEPOL 0.0 # a highly specialized CHARMM-related keyword: do not use unless you know exactly why") if (pmo.patching == 1) keyf.puts("FMCSC_CPATCHFILE ???? # partial charges are often context-dependent, e.g., in molecules with a mix of supported and unsupported residues") else keyf.puts("# FMCSC_CPATCHFILE ???? # partial charges often describe nonlocal (context-dependent) phenomena, e.g., terminal residues in AMBER") end if (pmo.ghosting == 1) keyf.puts("FMCSC_FEG_IPP 1.0 # scale factor for ghost version of inverse power potential; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_ATTLJ 1.0 # scale factor for ghost version of dispersion term; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_LJMODE 2 # cavity-forming potentials cannot just be scaled linearly: 2 is a buffered version\nFMCSC_FEG_LJRAD 0.5 # the buffer radius for the scaled potential in Angstrom\n# FMCSC_FEG_LJEXP ???? # this can be used to fine-tune the buffered potential\n# FMCSC_FEG_LJSCEXP ???? # another parameter to fine-tune the buffered potential\nFMCSC_FEG_POLAR 1.0# scale factor for the ghost version of the Coulomb potential\nFMCSC_FEG_CBMODE 1 # linear scaling is generally preferred for Coulomb\n# FMCSC_FEG_CBRAD ???? # would become relevant with buffered versions of the ghosted Coulomb potential\n# FMCSC_FEG_CBEXP ???? # ditto\n# FMCSC_FEG_CBSCEXP ???? # ditto") end need_lrel = 1 end when 3,4,5 pmo.stdnbs = 1 keyf.puts(ljs) if (pmo.patching == 1) keyf.puts("# FMCSC_LJPATCHFILE ???? # in some cases, it may prove useful to selectively patch LJ parameters (rather than biotypes)") end keyf.puts("FMCSC_SC_POLAR 1.0 # scale factor for Coulomb potential\nFMCSC_DIPREPORT 0 # can be used to produce VMD files for visualizing charge group information\n# FMCSC_AMIDEPOL 0.0 # a highly specialized CHARMM-related keyword: do not use unless you know exactly why") if (pmo.patching == 1) keyf.puts("FMCSC_CPATCHFILE ???? # partial charges are often context-dependent, e.g., in molecules with a mix of supported and unsupported residues") else keyf.puts("# FMCSC_CPATCHFILE ???? # partial charges often describe nonlocal (context-dependent) phenomena, e.g., terminal residues in AMBER") end if (pmo.ghosting == 1) keyf.puts("FMCSC_FEG_IPP 1.0 # scale factor for ghost version of inverse power potential; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_ATTLJ 1.0 # scale factor for ghost version of dispersion term; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_LJMODE 2 # cavity-forming potentials cannot just be scaled linearly: 2 is a buffered version\nFMCSC_FEG_LJRAD 0.5 # the buffer radius for the scaled potential in Angstrom\n# FMCSC_FEG_LJEXP ???? # this can be used to fine-tune the buffered potential\n# FMCSC_FEG_LJSCEXP ???? # another parameter to fine-tune the buffered potential\nFMCSC_FEG_POLAR 1.0# scale factor for the ghost version of the Coulomb potential\nFMCSC_FEG_CBMODE 1 # linear scaling is generally preferred for Coulomb\n# FMCSC_FEG_CBRAD ???? # would become relevant with buffered versions of the ghosted Coulomb potential\n# FMCSC_FEG_CBEXP ???? # ditto\n# FMCSC_FEG_CBSCEXP ???? # ditto") end need_lrel = 1 end if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) if (pmo.absinth == 1) keyf.puts("FMCSC_MOL2FOSLIBFILE ???? # you will need to supply a solvation group library file in conjunction with FMCSC_MOL2FILE to enable ABSINTH for small molecules") end keyf.puts("FMCSC_MOL2POLMODE 1 # for the screened molecules' input charges, mode 1 will force-regularize individual charge group total to net integers (input charges sometimes are not regularized at all, depending on their source)\nFMCSC_POLTOL 0.05 # a tolerance setting for determining groups of atoms comprising integer charges; the default of 0.0 will not work well for small molecule charges\n") else keyf.puts("# FMCSC_POLTOL 0.0 # a tolerance setting for determining groups of atoms comprising integer charges; may for example become relevant with patches or AMBER charges\n") end puts "Keywords related to standard nonbonded force field terms have been added.\n\n" when 2 keyf.puts(ljs) if (pmo.ghosting == 1) keyf.puts("FMCSC_FEG_IPP 1.0 # scale factor for ghost version of inverse power potential; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_ATTLJ 1.0 # scale factor for ghost version of dispersion term; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_LJMODE 2 # cavity-forming potentials cannot just be scaled linearly: 2 is a buffered version\nFMCSC_FEG_LJRAD 0.5 # the buffer radius for the scaled potential in Angstrom\n# FMCSC_FEG_LJEXP ???? # this can be used to fine-tune the buffered potential\n# FMCSC_FEG_LJSCEXP ???? # another parameter to fine-tune the buffered potential") end puts "Keywords related to Lennard-Jones potentials have been added.\n\n" when 3 keyf.puts("FMCSC_SC_WCA 1.0 # scale factor for Weeks-Chandler-Andersen potential\nFMCSC_ATT_WCA ???? # the WCA potential is used to interpolate the attractive term while keeping the minimum position constant: this is the scale for the dispersive well depth in units of epsilon_ij\nFMCSC_CUT_WCA ???? # the distance, in units of sigma_ij, where the potentials for atoms i,j becomes zero exactly\nFMCSC_SIGRULE ???? # please note that no supported force fields are parameterizaed with WCA in mind\nFMCSC_EPSRULE ???? # please note that no supported force fields are parameterizaed with WCA in mind\nFMCSC_HSSCALE 1.0 # can be used to scale sigma parameters globally") when 4 keyf.puts("FMCSC_SC_IPP 1.0 # scale factor for inverse power potential (excluded volume)\nFMCSC_SC_ATTLJ 0.0 # scale factor for dispersion potential (LJ)\nFMCSC_IPPEXP 12 # 12 is also the default (LJ)\nFMCSC_SIGRULE ???? # please look up correct value for your parameter file in documentation\nFMCSC_EPSRULE ???? # please look up correct value for your parameter file in documentation\nFMCSC_HSSCALE 1.0 # can be used to scale sigma parameters globally") if (pmo.ghosting == 1) keyf.puts("FMCSC_FEG_IPP 1.0 # scale factor for ghost version of inverse power potential; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_ATTLJ 0.0 # scale factor for ghost version of dispersion term; at 0.0, there is no contribution\nFMCSC_FEG_LJMODE 2 # cavity-forming potentials cannot just be scaled linearly: 2 is a buffered version\nFMCSC_FEG_LJRAD 0.5 # the buffer radius for the scaled potential in Angstrom\n# FMCSC_FEG_LJEXP ???? # this can be used to fine-tune the buffered potential\n# FMCSC_FEG_LJSCEXP ???? # another parameter to fine-tune the buffered potential") end when 5 keyf.puts("FMCSC_SC_WCA 1.0 # scale factor for Weeks-Chandler-Andersen potential\nFMCSC_ATT_WCA ???? # the WCA potential is used to interpolate the attractive term while keeping the minimum position constant: this is the scale for the dispersive well depth in units of epsilon_ij\nFMCSC_CUT_WCA ???? # the distance, in units of sigma_ij, where the potentials for atoms i,j becomes zero exactly\nFMCSC_SIGRULE ???? # please note that no supported force fields are parameterizaed with WCA in mind\nFMCSC_EPSRULE ???? # please note that no supported force fields are parameterizaed with WCA in mind\nFMCSC_HSSCALE 1.0 # can be used to scale sigma parameters globally") if (pmo.patching == 1) keyf.puts("# FMCSC_LJPATCHFILE ???? # in some cases, it may prove useful to selectively patch LJ parameters (rather than biotypes)") end keyf.puts("FMCSC_SC_POLAR 1.0 # scale factor for Coulomb potential\nFMCSC_DIPREPORT 0 # can be used to produce VMD files for visualizing charge group information\n# FMCSC_AMIDEPOL 0.0 # a highly specialized CHARMM-related keyword: do not use unless you know exactly why") if (pmo.patching == 1) keyf.puts("FMCSC_CPATCHFILE ???? # partial charges are often context-dependent, e.g., in molecules with a mix of supported and unsupported residues") else keyf.puts("# FMCSC_CPATCHFILE ???? # partial charges often describe nonlocal (context-dependent) phenomena, e.g., terminal residues in AMBER") end need_lrel = 1 if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) keyf.puts("FMCSC_MOL2POLMODE 1 # for the screened molecules' input charges, mode 1 will force-regularize individual charge group total to net integers (input charges sometimes are not regularized at all, depending on their source)\nFMCSC_POLTOL 0.05 # a tolerance setting for determining groups of atoms comprising integer charges; the default of 0.0 will not work well for small molecule charges\n") else keyf.puts("# FMCSC_POLTOL 0.0 # a tolerance setting for determining groups of atoms comprising integer charges; may for example become relevant with patches or AMBER charges\n") end when 6 keyf.puts("FMCSC_SC_IPP 1.0 # scale factor for inverse power potential (excluded volume)\nFMCSC_SC_ATTLJ 0.0 # scale factor for dispersion potential (LJ)\nFMCSC_IPPEXP 12 # 12 is also the default (LJ)\nFMCSC_SIGRULE ???? # please look up correct value for your parameter file in documentation\nFMCSC_EPSRULE ???? # please look up correct value for your parameter file in documentation\nFMCSC_HSSCALE 1.0 # can be used to scale sigma parameters globally") if (pmo.patching == 1) keyf.puts("# FMCSC_LJPATCHFILE ???? # in some cases, it may prove useful to selectively patch LJ parameters (rather than biotypes)") end keyf.puts("FMCSC_SC_POLAR 1.0 # scale factor for Coulomb potential\nFMCSC_DIPREPORT 0 # can be used to produce VMD files for visualizing charge group information\n# FMCSC_AMIDEPOL 0.0 # a highly specialized CHARMM-related keyword: do not use unless you know exactly why") if (pmo.patching == 1) keyf.puts("FMCSC_CPATCHFILE ???? # partial charges are often context-dependent, e.g., in molecules with a mix of supported and unsupported residues") else keyf.puts("# FMCSC_CPATCHFILE ???? # partial charges often describe nonlocal (context-dependent) phenomena, e.g., terminal residues in AMBER") end if (pmo.ghosting == 1) keyf.puts("FMCSC_FEG_IPP 1.0 # scale factor for ghost version of inverse power potential; at 1.0, the potential is equal to the standard one\nFMCSC_FEG_ATTLJ 0.0 # scale factor for ghost version of dispersion term; at 0.0, the potential is inactive\nFMCSC_FEG_LJMODE 2 # cavity-forming potentials cannot just be scaled linearly: 2 is a buffered version\nFMCSC_FEG_LJRAD 0.5 # the buffer radius for the scaled potential in Angstrom\n# FMCSC_FEG_LJEXP ???? # this can be used to fine-tune the buffered potential\n# FMCSC_FEG_LJSCEXP ???? # another parameter to fine-tune the buffered potential\nFMCSC_FEG_POLAR 1.0# scale factor for the ghost version of the Coulomb potential\nFMCSC_FEG_CBMODE 1 # linear scaling is generally preferred for Coulomb\n# FMCSC_FEG_CBRAD ???? # would become relevant with buffered versions of the ghosted Coulomb potential\n# FMCSC_FEG_CBEXP ???? # ditto\n# FMCSC_FEG_CBSCEXP ???? # ditto") end if (pmo.patching == 1) keyf.puts("FMCSC_CPATCHFILE ???? # partial charges are often context-dependent, e.g., in molecules with a mix of supported and unsupported residues") else keyf.puts("# FMCSC_CPATCHFILE ???? # partial charges often describe nonlocal (context-dependent) phenomena, e.g., terminal residues in AMBER") end need_lrel = 1 if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) keyf.puts("FMCSC_MOL2POLMODE 1 # for the screened molecules' input charges, mode 1 will force-regularize individual charge group total to net integers (input charges sometimes are not regularized at all, depending on their source)\nFMCSC_POLTOL 0.05 # a tolerance setting for determining groups of atoms comprising integer charges; the default of 0.0 will not work well for small molecule charges\n") else keyf.puts("# FMCSC_POLTOL 0.0 # a tolerance setting for determining groups of atoms comprising integer charges; may for example become relevant with patches or AMBER charges\n") end when 7 if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) keyf.puts("FMCSC_MOL2BONDMODE 3 # there are very few good reasons for not having this set to 3 in a small molecule screen\nFMCSC_PLANAR_TOLS 10.0 20.0 # the defaults are quite stringent") end keyf.puts("FMCSC_SC_BONDED_B 1.0 # scale factor for bond stretch potentials\nFMCSC_SC_BONDED_A 1.0 # scale factor for angle bending potentials\nFMCSC_SC_BONDED_T 1.0 # scale factor for dihedral angle potentials (usually cosine-based)\nFMCSC_SC_BONDED_I 1.0 # scale factor for improper torsional potentials (usually harmonic or quasi-harmonic)\nFMCSC_IMPROPER_CONV ???? # please look up the correct value for the force field in use\n# FMCSC_GUESS_BONDED 0 # turn this on to guess certain types of bonded potentials from input geometries (for unsupported residues in particular)\n# FMCSC_BONDREPORT 0 # can be turned on to obtain a complete summary of assigned/missing bonded potentials") if (pmo.patching == 1) keyf.puts("FMCSC_BPATCHFILE ???? # specific bonded terms can be patched using this; ") else keyf.puts("# FMCSC_BPATCHFILE ???? # specific bonded terms can be patched using this; for larger changes, biotype patches are better") end if (pmo.ghosting == 1) keyf.puts("FMCSC_FEG_BONDED_B 1.0 # scale factor for ghost version of bond stretch potentials\nFMCSC_FEG_BONDED_A 1.0 # scale factor for ghost version of angle bending potentials\nFMCSC_FEG_BONDED_T 1.0 # scale factor for ghost verion of dihedral angle potentials (usually cosine-based)\nFMCSC_FEG_BONDED_I 1.0 # scale factor for ghost version of improper torsional potentials (usually harmonic or quasi-harmonic") end puts "Keywords related to standard bonded potentials have been added.\n\n" when 8 keyf.puts("FMCSC_SC_BONDED_M 1.0 # scale factor for CMAP potentials (two consecutive dihedrals)\n# FMCSC_CMAPORDER ???? # this is only relevant if any CMAP potential uses the B-spline rather than the interpolating bicubic spline form") puts "Keywords related to CMAP potentials have been added.\n\n" when 9 keyf.puts("FMCSC_SC_TOR 1.0 # scale factor for dihedral angle-based bias potentials (parameter file-independent)\nFMCSC_TORFILE ???? # this is the required input file for this bias term\nFMCSC_TORREPORT 1 # reports are useful to enable when using custom file-based inputs") puts "Keywords related to torsional bias potentials have been added.\n\n" when 10 keyf.puts("FMCSC_SC_ZSEC 1.0 # scale factor for harmonic, torsion-based secondary structure bias potential\nFMCSC_ZS_FR_A ???? # target alpha content (0-1): only relevant if alpha restraint strength is nonzero\nFMCSC_ZS_FR_B ???? # target beta content (0-1): only relevant if beta restraint strength is nonzero\nFMCSC_ZS_FR_KA ???? # alpha restraint strength in kcal/mol\nFMCSC_ZS_FR_KB ???? # beta restraint strength in kcal/mol\n# The next 6 commented keywords will be listed below in uncommented form if you requested the corresponding analysis (FMCSC_SEGCALC)\n# FMCSC_ZS_POS_A ???? # you can move the alpha circle's center (phi/psi in deg)\n# FMCSC_ZS_POS_B ???? # you can move the beta circle's center (phi/psi in deg)\n# FMCSC_ZS_RAD_A ???? # you can change the alpha circle's radius (in deg)\n# FMCSC_ZS_RAD_B ???? # you can change the beta circle's radius (in deg)\n# FMCSC_ZS_STP_A ???? # you can change the rate of decay outside of the alpha circle (in 1/deg^2)\n# FMCSC_ZS_STP_B ???? # you can change the rate of decay outside of the beta circle (in 1/deg^2)\n") puts "Keywords related to torsional secondary structure restraints have been added.\n\n" when 11 keyf.puts("FMCSC_SC_DSSP 1.0 # scale factor for harmonic, hydrogen bond-based secondary structure bias potential\nFMCSC_DSSP_HSC ???? # target alpha content (0-1): only relevant if alpha restraint strength is nonzero\nFMCSC_DSSP_ESC ???? # target beta ('extended') content (0-1): only relevant if beta restraint strength is nonzero\nFMCSC_DSSP_HSC_K ???? # alpha restraint strength in kcal/mol\nFMCSC_DSSP_ESC_K ???? # beta restraint strength in kcal/mol\nFMCSC_DSSP_MODE ???? # since the raw content variables are discrete, the restraint is applied to a continuous transform, the specific form of which is set by this keyword\n# the next 5 keywords will be listed below (3 in uncommented form) if you are using the DSSP analysis feature\n# FMCSC_DSSP_EXP ???? # see FMCSC_DSSP_MODE\n# FMCSC_DSSP_GOODHB ???? # you can change the default (in kcal/mol) for what is a good hydrogen bond\n# FMCSC_DSSP_MAXHB ???? # you can change the threshold (in kcal/mol) for what to consider as a hydrogen bond\n# FMCSC_DSSP_MINHB ???? # you can change the default (in kcal/mol) for the minimum (best) possible energy of any hydrogen bond\n# FMCSC_DSSP_CUT ???? # you can change the default cutoff (in Angstrom) used for searching H-bonds") puts "Keywords related to DSSP secondary structure restraints have been added.\n\n" when 12 keyf.puts("FMCSC_SC_POLY 1.0 # scale factor for polymeric bias potentials to act on the sizes and shapes of molecules\nFMCSC_POLYFILE ???? # this is the required input file for this bias term\nFMCSC_POLYREPORT 1 # reports are useful to enable when using custom file-based inputs") puts "Keywords related to polymeric restraints have been added.\n\n" when 13 keyf.puts("# Note that tabulated interactions (TABUL) are considered long-range nonbonded terms but have no long-range correction scheme in place.\n# This means that cutoffs may incur significant artifacts unless the potentials are chosen appropriately.\nFMCSC_SC_TABUL 1.0 # scale factor for arbitrary interatomic potentials acting on pairs of atoms\nFMCSC_TABCODEFILE ???? # this is one of the required input files for this bias term: it determines the map of potentials to atom pairs\nFMCSC_TABPOTFILE ???? # this is one of the required input files for this bias term: it defines the potentials as a function of distance (kcal/mol and Angstrom)\n# FMCSC_TABTANGFILE ???? # you can choose to supply potential gradients by file as well rather than taking finite differences (kcal/mol/Angstrom)\n# FMCSC_TABITIGHT ???? # you could manipulate the finite difference inference in special cases\n# FMCSC_TABIBIAS ???? # you could manipulate the finite difference inference in special cases\nFMCSC_TABREPORT 1 # reports are useful to enable when using custom file-based inputs") puts "Keywords related to tabulated interatomic potentials have been added.\n\n" when 14 keyf.puts("FMCSC_SC_DREST 1.0 # scale factor for (half-)harmonic distance restraints to act on pairs of atoms\nFMCSC_DRESTFILE ???? # this is the required input file for this bias term\nFMCSC_DRESTREPORT 1 # reports are useful to enable when using custom file-based inputs") puts "Keywords related to custom distance restraints have been added.\n\n" when 15 keyf.puts("FMCSC_SC_OSMO 1.0 # scale factor for half-harmonic confinement (compartmentalization) potentials to act on molecules or molecule types\nFMCSC_OSMOFILE ???? # this is the required input file for this bias term\nFMCSC_OSMO_MODE ???? # this keyword sets (indirectly) the number of compartments to be created\nFMCSC_OSMOREPORT 1 # reports are useful to enable when using custom file-based inputs") puts "Keywords related to compartmentalization potentials have been added.\n\n" when 16 keyf.puts("FMCSC_SC_EMICRO 1.0 # scale factor for lattice-based spatial density restraints\nFMCSC_EMMAPFILE ???? # this is the required input file for this bias term; it contains the 3D lattice with the target density in a NetCDF format\nFMCSC_EMMODE 1 # changing this means you are switching to an average rather rhan instantanoues density restraint (experimental)\n# FMCSC_EMIWEIGHT ???? # would be needed for working with restraints on time averages\n# FMCSC_EMREDUCE ???? ???? ???? # it is possible to reduce the resolution of the input map (target xyz cell sizes in Angstrom)\n# FMCSC_EMBACKGROUND ???? # you can use this to override the automatically determined background level in the input map\n# FMCSC_EMFLATTEN ???? # this keyword can be used to cut density spikes in the input to this value\n# FMCSC_EMTRUNCATE ???? # this keyword can be used to raise density troughs in the input map to this value\nFMCSC_EMTHRESHOLD ???? # this sets the density level above which the true signal should be started\nFMCSC_EMTOTMASS ???? # this sets the total mass that should explain the true signal in g/mol\nFMCSC_EMHEURISTIC 3 # for asymmetric geometries, the simpler heuristics (choices 1 or 2) may be more efficient\n# FMCSC_EMPROPERTY 4 # the default is to use mass, choice 1\n# FMCSC_EMPROPERTYFILE ???? # if EMPROPERTY is 4, you can supply a custom property value per atom by file") puts "Keywords related to lattice-based spatial density restraints have been added.\n\n" end if ((choice <= 4) && (choice >= 2) && (pmo.patching == 1)) keyf.puts("# FMCSC_LJPATCHFILE ???? # in some cases, it may prove useful to selectively patch LJ parameters (rather than biotypes)") end end if (terms_selected.include?(13)) pmo.stdnbs = 0 end if (need_lrel == 1) nest2_lrel(keyf,pmo) end if (pmo.ghosting == 1) if ( (!terms_selected.include?(2)) && (!terms_selected.include?(4)) && (!terms_selected.include?(6)) && (!terms_selected.include?(7)) && (!terms_selected.include?(1)) ) puts "Please note: you indicated to use ghosting but did not select any energy term that would support this functionality." keyf.puts("FMCSC_GHOST 0 # no ghosting of (primarily nonbonded) interactions will occur; above FEG-related keywords are irrelevant") end elsif ((terms_selected.length > 0) && (pmo.dockrun != 1) && (pmo.docksrun != 1)) puts "The final Hamiltonian-related question is whether you want to use energy landscape sculpting, essentially the technique originally made famous as 'Accelerated Molecular Dynamics' (y/n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("# Energy landscape sculpting changes the potential energy landscape at the level of individual terms.\n# The sculpting consists of shaving barriers and filling basins, and eligible terms are those enabled.\n# Limitations exist regarding nonbonded terms (not separable in dynamics) and total energy (mutually exclusive with all other terms)\nFMCSC_SCULPT ???? # possible choices generally correspond to column numbers of ENERGY.dat; plus '15' for total nonbonded terms\nFMCSC_ELS_PRINT_WEIGHTS 10 # interval for reporting snapshot weights in ELS sense: allows diagnosis of relevance of ensemble\nFMCSC_ELS_FILLS ???? # as many fill values (kcal/mol) as terms for FMCSC_SCULPT\nFMCSC_ELS_ALPHA_F ???? # buffer values (positive, kcal/mol) for basin fills (larger means smoother + weaker)\nFMCSC_ELS_SHAVES ???? # as many shave values (kcal/mol) as terms for FMCSC_SCULPT\nFMCSC_ELS_ALPHA_S ???? # buffer values (positive, kcal/mol) for barrier shaves (larger means smoother + weaker)") else keyf.puts("# FMCSC_SCULPT 0 # it is safest that this keyword is simply not present if you do not want to use energy landscape sculpting\n") end end keyf.puts("#\n# FMCSC_CMAPDIR ???? # usually the CAMPARI data/ subdirectory: you will need this keyword for any parameter file that prescribes CMAP terms irrespective of whether you use CMAP or not in the actual run\n") end def nest1_analysis(keyf,pmo,isa) keyf.puts("#\n# Snapshot-by-snapshot analysis choices") if (isa == 1) puts "One of the main purposes of CAMPARI's trajectory analysis mode is to perform snapshot-by-snapshot analyses. " base1 = 1 base2 = 1 else puts "CAMPARI can perform snapshot-by-snapshot analyses on-the-fly for a running simulation. This is not only simpler than post-processing but can also be useful by utilizing more data than one is willing store permanently." base1 = 50 base2 = 500 end if ((pmo.pmode == 1) || (pmo.pmode == 3)) keyf.puts("# Be aware that a proper load balance when using threads can be difficult to achieve for these tasks (see documentation).") end terms_unselected = Array.new(12){ |i| i+1 } choice = 0 while (choice != 13) puts "\nWhat snapshot-by-snapshot analyses do you want to perform?" terms_unselected.each do |i| printf("(%2d): ",i) case i when 1 puts "Cheap generic polymer properties (size and shape histograms. radial density profiles, turn angle averages along polymer backbones and their correlation functions)" when 2 puts "Expensive generic polymeric properties (internal scaling behavior of polymers, estimates of hydrodynamic radii, scattering profiles as obtained in SAXS experiments)" when 3 puts "Contacts analysis (contact maps, contact number histograms, contact-based assessment of the cluster formation of molecules)" when 4 puts "Polypeptide backbone dihedral statistics (Ramachandran maps, vicinal 1HA-1HN NMR coupling constants)" when 5 puts "Dihedral angle-based secondary structure analysis for polypeptides (histograms and segment length distributions as used for example in helix-coil models)" when 6 puts "Hydrogen bond-based secondary structure analysis for polypeptides (DSSP, includes histograms, segment length distributions and the option to obtain the full instantaneous DSSP output)" when 7 puts "Spatial density analysis (distribution of molecular properties (like mass) in absolute 3D space (on a lattice))" when 8 puts "Pair correlation (distance statistics) analysis (can be used to also obtain traces of instantaneous distances, e.g., to compute FRET efficiencies)" when 9 puts "General statistics on molecular coordinates (topology-derived angles, dihedral angles, etc.)" when 10 puts "Solvent-accessible volume analysis (solvent accessibility is calculated as in the ABSINTH implicit solvation model)" when 11 puts "Dipole analysis (mean dipole moments at various resolution levels)" when 12 puts "Custom analysis through the Python interface (you can do whatever you want here, see Tutorial 18)" end end puts "(13): I do not want to perform any of these analyses or I have selected all the ones I want." choice = curated_int_answer(gets.chomp,1,13) if (choice < 13) terms_unselected.delete(choice) end case choice when 1 keyf.print("FMCSC_POLCALC ",base1," # computes polymeric properties that depend linearly on polymer size (cheap)\nFMCSC_RGBINSIZE 0.2 # bin size in Angstrom for RGHIST.dat, RETEHIST.dat, and DENSPROF.dat\nFMCSC_POLRGBINS 1000 # this is the default number of bins\n") puts "Keywords related to computationally inexpensive polymer analyses have been added.\n\n" when 2 keyf.puts("FMCSC_RHCALC ???? # these analyses scale quadratically with polymer size and thus become very expensive for chains with 100s of residues\nFMCSC_SCATTERCALC ???? # this interval is relative to FMCSC_RHCALC\nFMCSC_SCATTERRES ???? # only relevant if FMCSC_SCATTERCALC*FMCSC_RHCALC < FMCSC_NRSTEPS\nFMCSC_SCATTERVECS ???? # only relevant if FMCSC_SCATTERCALC*FMCSC_RHCALC < FMCSC_NRSTEPS") puts "Keywords related to computationally expensive polymer analyses have been added.\n\n" when 3 keyf.puts("FMCSC_CONTACTCALC ???? # scales quadratically with the number of solute residues\nFMCSC_CONTACTOFF 0 # there is no obvious benefit of using numbers larger than 0 here\nFMCSC_CONTACTMIN ???? # atom-atom based contact criterion (any pair)\nFMCSC_CONTACTCOM ???? # center of mass-based contact criterion\nFMCSC_CLUSTERCALC ???? # interval is relative to FMCSC_CONTACTCALC, implied contact def. relies on FMCSC_CONTACTMIN") puts "Keywords related to contact and solution cluster analyses have been added.\n\n" when 4 keyf.print("FMCSC_ANGCALC ",base1," # computes Ramachandran maps and NMR coupling constants (cheap)\nFMCSC_ANGRES 6.0 # bin size in degrees for RAMACHANDRAN.dat, RAMARES_*dat, and RAMAMOL*.dat\nFMCSC_RAMARES 0 # can be used to request maps for a few specific residues (number is first, followed by indices)\nFMCSC_RAMAMOL 0 # can be used to request maps averaged for a few specific molecules (number is first, followed by indices)\n") puts "Keywords related to polypeptide backbone analyses of dihedral angles have been added.\n\n" when 5 keyf.print("FMCSC_SEGCALC ",base1," # computes segment statistics and secondary structure histograms (cheap)\nFMCSC_BBSEGFILE ????/data/bbseg2.dat # the file used for the Ramchandran discretization underlying segments statistics in BB_SEGMENTS*\nFMCSC_ZSEC_POS_A -60. 50. # suitable as central position for alpha (deg)\nFMCSC_ZSEC_POS_B -155. 160. # suitable as central position for antiparallel beta (deg)\nFMCSC_ZSEC_RAD_A 35. # radius of alpha-circle (deg)\nFMCSC_ZSEC_RAD_B 35. # radius of beta circle (deg)\nFMCSC_ZSEC_STP_A 0.002 # default\nFMCSC_ZSEC_STP_B 0.002 # default\n") puts "Keywords related to polypeptide secondary structure analyses based on dihedral angles have been added.\n\n" when 6 keyf.print("FMCSC_DSSPCALC ",base2," # DSSP statistics have to be calculated globally and use cutoffs to tame the quadratic complexity with the number of H-bond sites\n# FMCSC_DSSP_MODE ???? # there are different ways to get weighted DSSP scores; this keyword will be listed above already if the DSSP restraint potential was selected\n# FMCSC_DSSP_EXP ???? # see FMCSC_DSSP_MODE\nFMCSC_DSSP_GOODHB -2.5 # you can set the energy in kcal/mol of a good hydrogen bond (does not affect assignment per se, matter only for weighted scores)\nFMCSC_DSSP_MAXHB -0.5 # you can change the threshold (in kcal/mol) for what to consider as a hydrogen bond; this is the default and the DSSP reference value (affects everything)\nFMCSC_DSSP_MINHB -5.0 # you can set the minimum (best) possible energy in kcal/mol of any hydrogen bond (does not affect assignment per se, matter only for weighted scores)\n# FMCSC_DSSP_CUT ???? # you could change the default cutoff (in Angstrom) used for searching H-bonds (should not be needed)\n") if ((pmo.mpi_mode != 1) || ((pmo.pmode != 2) && (pmo.pnode != 3))) puts "Do you require instantaneous DSSP output (output files can become very large), which is the only way of knowing the details of each instantaneous hydrogen bond pattern (y/n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_INSTDSSP 1 # 1 enables writing of every analyzed, instantaneous DSSP pattern to DSSP_RUNNING.dat") else keyf.puts("FMCSC_INSTDSSP 0 # 0 disables writing of instantaneous DSSP (to DSSP_RUNNING.dat)") end end puts "Keywords related to DSSP analysis have been added.\n\n" when 7 keyf.puts("FMCSC_EMCALC ???? # the cost of spatial density analysis depends strongly on the number of relevant particles in the system and the number of lattice cells; it is free if the restraint potential is on\nFMCSC_EMDELTAS ???? ???? ???? # lattice spacing (x,y,z)\nFMCSC_EMBSPLINE ???? # the B-spline interpolation order: the number of lattice cells each particle 'touches' grows cubically with this order\nFMCSC_EMBGDENSITY 0.0 # the assumed background level should be different when in an implicit solvent\nFMCSC_EMBUFFER 1.2 # irrelevant only in full 3D periodic boundary conditions") puts "Keywords related to spatial density analysis have been added.\n\n" when 8 puts "This covers different functionalities with solute-solute (center-of-mass) distance monitoring always enabled (by pair of types). Do you also want to monitor specific atom-atom distances (y/n)?" keyf.puts("# The default component of the pair correlation facility can be expensive.\n# The cost can be regulated using FMCSC_ANGRPFILE to tag molecule types as solvent (but this extends then to other analyses)\nFMCSC_PCCALC ???? # the cost of these analyses is linear in the number of monitored distances(!) (no cutoffs or simplifications apply)\nFMCSC_PCBINSIZE 0.2 # the bin size for the histograms in Angstrom\nFMCSC_DO_AMIDEPC 0 # 0 disables intramolecular analysis of distances between sidechain amides in peptides") if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_PCCODEFILE ???? # required custom input file for selecting atom-atom distance pairs\nFMCSC_GPCREPORT 1 # the report writes the file 'GENERAL_PC.idx'") puts "Do you require instantaneous distance traces, e.g., for cross-correlation analyses (y/n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_INSTGPC ???? # relative (to FMCSC_PCCALC) frequency for writing instantaneous atom-atom distances to GENERAL_DIS.dat") else keyf.puts("FMCSC_INSTGPC 0 # 0 disables writing of any instantaneous atom-atom distances (to GENERAL_DIS.dat)") end end puts "Keywords related to pair correlation (distance distribution) analysis have been added.\n\n" when 9 keyf.print("FMCSC_INTCALC ",base1," # increments Z-matrix coordinate histograms (cheap)\nFMCSC_WHICHINT 0 0 0 1 # this enables only dihedral angles (INTHISTS_DI.dat); order is bonds, angles, impropers, dihedrals\nFMCSC_TOROUT 2000000000 # change this to a small negative or positive value to get instantanteous data for modifiable/native dihedral angles to FYC.dat\nFMCSC_COVCALC 2000000000 # use would be redundant\n") puts "Keywords related to internal coordinate analysis have been added.\n\n" when 10 if (pmo.absinth == 1) keyf.print("# Note that you cannot use independent parameters for the SAV estimation in the DMFI and the analysis functionality\nFMCSC_SAVCALC ",base1," # the cost is negligible when ABSINTH is turned on\n") if (pmo.patching == 0) keyf.puts("# FMCSC_ASRPATCHFILE ???? # could be used to override atomic volume reduction factors per atom\n# FMCSC_SAVPATCHFILE ???? # could be used to override maximum solvent accessible volume fractions per atom\n# FMCSC_RPATCHFILE ???? # could be used to override atomic radii per atom") end else keyf.puts("# Note that you cannot use independent parameters for the SAV estimation in the DMFI and the analysis functionality\nFMCSC_SAVCALC ???? # this would be handled more efficiently if the ABSINTH model were turned on (FMCSC_SC_IMPSOLV)\nFMCSC_SAVPROBE ???? # sets the 'probe' radius (half the thickness of the volume shell)") keyf.puts("# FMCSC_ASRPATCHFILE ???? # could be used to override atomic volume reduction factors per atom\n# FMCSC_SAVPATCHFILE ???? # could be used to override maximum solvent accessible volume fractions per atom\n# FMCSC_RPATCHFILE ???? # could be used to override atomic radii per atom") end puts "Do you require instantaneous SAV output, which is the only way of retrieving instantaneous, atomic SAVs along a trajectory (y/n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_INSTSAV ???? # relative (to FMCSC_SAVCALC) frequency for writing of instantaneous SAVs globally and for selected atoms to SAV.dat\nFMCSC_SAVATOMFILE ???? # instantaneous SAVs are useful primarily for individual atoms of interest (custom input file)") else keyf.puts("FMCSC_INSTSAV 0 # 0 disables writing of any instantaneous SAVs (to SAV.dat)") end puts "Keywords related to solvent-accessible volume analysis have been added.\n\n" when 11 keyf.print("FMCSC_DIPCALC ",base1," # compute dipole statistics (cheap); this requires the Coulomb potential to be turned on\n# FMCSC_SC_POLAR 1.0 # may need to uncomment\n") puts "Keywords related to dipole analysis have been added (this requires turning on an energy term if it is not already turned on).\n\n" when 12 keyf.print("# The Python interface relies on a module file campana.py that must be findable and importable for ForPy (requires NumPy and\n# SciPy, but can import other module if you need them\n# You should use the same Python version that was used to link CAMPARI to and you should avoid\n# Python environments because low-level system library versions must be consistent (otherwise cryptic errors are produced)\nFMCSC_PYCALC ",base1," # call the Python interface: every time the function campana_step in campana.py will be called\nFMCSC_PYDIR . # normally, you make a copy of campana.py and put in the local working directory, which this should point to\n") end end end def nest1_insta(keyf,pmo,isa) puts "\nCAMPARI produces limited instantaneous output during a running simulation, most prominently trajectory data. These files can in general get rather large. Some instantaneous writes are covered elsewhere (keywords FMCSC_INSTGPC, FMCSC_INSTSAV, FMCSC_INSTDSSP, and FMCSC_TOROUT)." puts "Regarding trajectory output. What do you require?\n(1): No output\n(2): A pdb-trajectory (large files!, MODEL/ENDMDL syntax)\n(3): A dcd-trajectory (binary, sequential access only, CHARMM/NAMD)\n(4): An xtc-trajectory (binary, sequential access only, smallest files, GROMACS)\n(5): A NetCDF trajectory (binary, random access, AMBER)\n(6): Writing to a SQL database (best meta-information, random access, database must be created independently)\n" choice = curated_int_answer(gets.chomp,1,6) if (choice > 1) puts "And what parts of the system should be written to these files (including *START.pdb and *END.pdb)?\n(1): All of it\n(2): Everything labeled as solutes (see FMCSC_ANGRPFILE for control in this regard)\n(3): A custom selection of atoms (by index)\n(4): Solutes and a proximity-based, automatic selection of non-solute molecules\n" choice2 = curated_int_answer(gets.chomp,1,4) end keyf.puts("#\n# Fundamental analysis keywords\n# FMCSC_ANGRPFILE ???? # you will need to uncomment and put a valid file if you want to split identical molecules into different groups for analysis or if solvent molecules need to be included for some tasks\nFMCSC_DISABLE_ANALYSIS 1 # for clarity, let us manually disable all analysis to begin with (this is a meta-keyword that must appear after FMCSC_NRSTEPS!)") keyf.puts("#\n# Instantaneous output settings\n# FMCSC_PDB_OUTPUTSTRING \"a6,i5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3\" # this allows you to customize the format of all pdb files written\n# FMCSC_PDB_NUCMODE 2 # choosing the nondefault '2' asks CAMPARI to write nucleic acids such that the grouping of atoms is in pdb convention\n# FMCSC_PDB_W_CONV ???? # this can be used to improve the compatibility of CAMPARI-generated pdbs with other software") case choice when 1 keyf.puts("# FMCSC_XYZOUT ???? # trajectory output is disabled\n# FMCSC_XYZPDB ???? # this would set the file format") when 2 keyf.puts("FMCSC_XYZOUT ???? # the step interval to append the trajectories; small values can produce large files\nFMCSC_XYZPDB 2 # 2 is the default and chooses the pdb convention as the file format\nFMCSC_XYZMODE 2 # 2 is the default and relevant only for pdb-files") when 3 keyf.puts("FMCSC_XYZOUT ???? # the step interval to append the trajectories; small values can produce large files\nFMCSC_XYZPDB 3 # 3 chooses the CHARMM/NAMD dcd convention as the file format") when 4 keyf.puts("FMCSC_XYZOUT ???? # the step interval to append the trajectories; small values can produce large files\nFMCSC_XYZPDB 4 # 4 chooses the GROMACS-like xtc convention as the file format\n# FMCSC_XTCPREC ???? # an auxiliary keyword controlling xtc-precision (default is 1000., which is ~3 digits )") when 5 keyf.puts("FMCSC_XYZOUT ???? # the step interval to append the trajectories; small values can produce large files\nFMCSC_XYZPDB 5 # 5 chooses the AMBER NetCDF convention as the file format") when 6 keyf.puts("FMCSC_XYZOUT ???? # the step interval to append the trajectories; small values can produce large files\nFMCSC_XYZPDB 6 # 6 chooses the PostgreSQL database convention as the file format\nFMCSC_PSQL_W_TRAJKEY ???? # name of the data set in the database: this is a key that identifies specific systems (define or reuse if possible/desirable)\nFMCSC_PSQL_W_DATABASE ???? # the name of the database (write permissions required)\nFMCSC_PSQL_W_DBHOST ???? # the database host's (IP) address, default is 127.0.0.1 (local)\n# FMCSC_PSQL_W_DBPORT 5432 # change only if your database administrators changed the default port\nFMCSC_PSQL_SIMSUM ???? # file with text providing meta-information: if not provided, it will be the summary portion of the CAMPARI run\n# FMCSC_PSQL_W_SOFTWARE ???? # adjust if you are depositing a simulation created by different software\n# FMCSC_PSQL_W_SIMUOFF ???? # relevant when appending existing database entries") end if (choice > 1) if (choice2 == 1) keyf.puts("FMCSC_XYZ_SOLVENT 1 # 1 makes sure that molecules tagged as solvent (FMCSC_ANGRPFILE) will be included") elsif (choice2 == 2) if (isa == 0) keyf.puts("# Keep in mind that not including solvent coordinates in simulation output means that they are lost for good.") end keyf.puts("FMCSC_XYZ_SOLVENT 0 # 0 makes sure that molecules tagged as solvent (FMCSC_ANGRPFILE) are excluded") elsif (choice2 ==3) if (isa == 0) keyf.puts("# Keep in mind that all excluded coordinates are lost for good.") end keyf.puts("FMCSC_TRAJIDXFILE ???? # a required custom index file containing the atoms to be included in trajectory output") elsif (choice2 ==3) if (isa == 0) keyf.puts("# Keep in mind that all excluded coordinates are lost for good.") end keyf.puts("FMCSC_XYZ_PRUNE ???? # a required list of molecule types and numbers to retain\n# FMCSC_TRAJIDXFILE ???? # an optional reference set to define the criterion of proximity for retaining the N closest molecules of all selected types\nFMCSC_XYZ_SOLVENT 1 # 0 is possible but confusing in this mode") end end pmo.talign = 0 if (isa == 1) puts "In trajectory analysis mode, it is possible to align coordinates before proceeding with analysis. This will affect analysis features relying on absolute coordinates but can also influence relative coordinates. What do you require?\n(0): No alignment\n(1): Alignment to a fixed structure\n(2): Alignment to the preceding structure\n(3): Alignment to a fixed structre and computation of distances according to a different set of atoms\n(4): Alignment to the preceding structure and computation of distances according to a different set of atoms" pmo.talign = curated_int_answer(gets.chomp,0,3) if (pmo.talign > 0) keyf.puts("FMCSC_ALIGNCALC -1 # -1 enables alignment at every step and asks CAMPARI to print instantaneous RMSDs (to RMSD.dat)\nFMCSC_ALIGNFILE ???? # a required custom input file defining the set of atoms (by index) to align on\n# FMCSC_XYZ_REFMOL ???? # this keyword can be used to define a reference molecule for alignment and for placing the 'best' periodic image of every molecule on output") end if ((pmo.talign == 1) || (pmo.talign == 3)) if (pmo.havetmpl == 1) keyf.puts("# The use of FMCSC_PDB_TEMPLATE as the reference structure for alignment means that the corresponding file may have to fulfill up to 3 functions simultaneously.\n# FMCSC_PDB_TEMPLATE ???? # see above: a single file for all functions must suffice for FMCSC_PDB_TEMPLATE") else keyf.puts("FMCSC_PDB_TEMPLATE ???? # the PDB file containing the reference structure for alignment (must be complete)") end elsif ((pmo.talign == 3) || (pmo.talign == 4)) keyf.puts("FMCSC_CFILE ???? # supplying this custom input file allows the alignment set to be split from the distance set") end end if (choice > 1) if (pmo.talign == 0) keyf.puts("# FMCSC_XYZ_REFMOL ???? # this keyword can be used to define a reference molecule for placing the 'best' periodic image of every molecule on output") end if (pmo.cocho != 2) keyf.puts("# If you want to force all atoms to always reside in the central unit cell, you'll need to change FMCSC_XYZ_FORCEBOX above (if there) or uncomment below?") keyf.puts("# FMCSC_XYZ_FORCEBOX ???? # 2 would force atoms into unit cell while assuming whole molecules on any possible input; 3 would also assume broken on input") end end if (isa == 0) if (pmo.mcflag == 1) keyf.print("FMCSC_ACCOUT ",[1,(pmo.nsteps/100.0).to_i].max," # net MC move acceptance counts written to ACCEPTANCE.dat\n") end if ((pmo.simtype != 1) && (pmo.simtype != 6)) keyf.print("FMCSC_ENSOUT ",[100,[1,(pmo.nsteps/10000.0).to_i].max].min," # instantaneous ensemble information written to ENSEMBLE.dat\n") end if (pmo.simtype == 1) keyf.print("FMCSC_CHECKFREQ ",[100,(pmo.nsteps/1000.0).to_i].max," # report system ensemble variables infrequently in log-output (pure MC)\n") elsif (pmo.simtype == 6) keyf.print("FMCSC_CHECKFREQ ",[1,(pmo.nsteps/100000.0).to_i].max," # report system ensemble variables frequently in log-output (minimization)\n") else keyf.print("FMCSC_CHECKFREQ ",[10,(pmo.nsteps/10000.0).to_i].max," # report system ensemble variables in log-output\n") end end if (pmo.emptyh == 1) if (pmo.cocho != 1) keyf.print("# FMCSC_ENOUT ???? # with these settings, should write only the boundary potential\n") end else if (isa == 1) keyf.print("FMCSC_ENOUT ",1," # instantaneous energy terms written to ENERGY.dat\n") else keyf.print("FMCSC_ENOUT ",[1000,[1,(pmo.nsteps/10000.0).to_i].max].min," # instantaneous energy terms written to ENERGY.dat\n") end end if (pmo.mpi_mode == 2) puts "Since this is a multi-copy MPI run with possibly different Hamiltonians for different copies, there is the potential to compute energy-based overlap metrics across replicas. What do you require?\n(0): Nothing\n(1): Compute only average overlap metrics across neighboring replicas\n(2): Compute only average overlap metric across all replicas\n(3): Compute average and print instantaneous overlap metric across neighboring replicas\n(4): Compute average and print instantaneous overlap metric across all replicas\n" choice3 = curated_int_answer(gets.chomp,0,4) if (choice3 == 0) keyf.puts("# FMCSC_REOLCALC ???? # this would allow you to enable overlap metrics across replicas") elsif ((choice3 == 1) || (choice3 == 3)) keyf.puts("FMCSC_REOLCALC ???? # these global energy evaluations can be expensive (output to N_*_OVERLAP.dat)\nFMCSC_REOLALL 0 # 0 chooses neighboring replicas only") elsif ((choice3 == 2) || (choice3 == 4)) keyf.puts("FMCSC_REOLCALC ???? # the cost of the global energy evaluations is now linear in the number of replicas (output to N_*_OVERLAP.dat)\nFMCSC_REOLALL 1 # 1 chooses all replicas (most overlap measures will often be ~0.0)") end if ((choice3 == 3) || (choice3 == 4)) keyf.puts("FMCSC_REOLINST ???? # interval *relative* to FMCSC_REOLCALC how often to print instantaneous values, e.g., required for Bennett acceptance ratio method") end if (isa == 1) puts "If such a run is an analysis run, there is the option to read back the trace of a replica exchange calculation to 'demix' the trajectories. Does this apply to you (y/n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("# In general, the reading of MPI trace files refers to absolute simulation steps (first column) of the original simulation.\nFMCSC_TRACEFILE ???? # the name of the trace file used to demix REX trajectories\nFMCSC_RE_TRAJOUT ???? # this must correspond to the original trajectory output interval\nFMCSC_RE_TRAJSKIP ???? # this must correspond to the number of original steps for which trajectory output was suppressed at the beginning\nFMCSC_RE_TRAJTOTAL ???? # this does not refer to steps: it must the number of stored snapshots in each trajectory file") end end end end def nest2_mpiism(keyf,pmo) puts "\nIn MPI mode, CAMPARI generates multiple copies of a simulation or analysis system and task(s). Which of the two fundamental execution modes do you want?\n(1): All copies are identical, and information is pooled (used for example to increase sampling by independent copies, to speed up analyses of large data sets, or in PIGS).\n(2): All copies retain independent information (used for example in replica exchange simulations and their analyses and in distributed small molecule screens)." pmo.mpi_mode = curated_int_answer(gets.chomp,1,2) case pmo.mpi_mode when 1 keyf.puts("# MPI settings\nFMCSC_MPIAVG 1 # 1 enables MPIAVG parallel mode, which implies that copies pool their data\n# FMCSC_REPLICAS ???? # the number of copies is taken from the MPI environment if there is a mismatch anyway -> redundant keyword in this mode\n# FMCSC_RETRACE 1 # there is usually no reason to disable the writing of any possible trace file\n") when 2 keyf.puts("# MPI settings\nFMCSC_REMC 1 # 1 enables REX-style parallel mode, which implies that copies keep their data separated\nFMCSC_REPLICAS ???? # this is a required keyword to indicate the number of copies (MPI processes)\n") end keyf.puts("# FMCSC_PDB_MPIMANY 1 # turn this on if you have a simulation task with PDB input and want to use separate input files per replica\n# FMCSC_MPICOLLS 0 # if MPI performance is unsatisfactory, you could try switching to 1, i.e., internal functions (this is unlikely to help)\n# FMCSC_MPIGRANULESZ ???? # for the customization of CAMPARI-internal collective MPI routines") end def nest1_strucinit(keyf,pmo) if ((pmo.dockrun == 1) || (pmo.docksrun == 1)) puts "For starting a small molecule screen, CAMPARI should generally read the fixed part of the system from a file. What scenario best applies to your case?\n(1): The entire system (all atoms) should be read from PDB, and all present residues except the screened small molecules are CAMPARI-supported.\n(2): The entire system including additional unsupported residues (but excepting the screened small molecules) should be read from a single, complete PDB.\n(3): There is no fixed part of the system.\n(4): Something else." choice = curated_int_answer(gets.chomp,1,4) keyf.puts("#\n# Settings for how to generate the conformation the system initially resides in\n# Note that not all missing atoms in input PDB files involve relevant degrees of freedom.") case choice when 1,2 keyf.puts("FMCSC_PDBFILE ???? # if this input file is missing or cannot be read, CAMPARI will default to randomization!\nFMCSC_PDB_READMODE 3 # 3 is needed for CAMPARI to ignore the missing placeholder (should find complete input otherwise)\nFMCSC_PDB_HMODE 1 # no need to rebuild\nFMCSC_RANDOMIZE 0 # this option eliminates randomization wherever feasible") when 3 keyf.puts("FMCSC_PDBFILE ???? # if there is nothing else in the system, CAMPARI needs to find a mock PDB file containing one atom with the right residue name (otherwise it exits)") when 4 keyf.puts("FMCSC_PDBFILE ???? # if this input file is missing or cannot be read, CAMPARI will default to randomization!\nFMCSC_PDB_TEMPLATE ???? # CAMPARI expects to find sufficient information for unsupported residues here (see docu)\nFMCSC_PDB_HMODE 1 # try changing to 2 if problems arise with hydrogen positions after read-in\nFMCSC_PDB_READMODE 3 # 3 is default and instructs CAMPARI to use whatever possible from a PDB (also see tool 'convert_SEQRES_toseq.sh')") end if (choice != 3) keyf.puts("# FMCSC_PDB_INPUTSTRING \"a6,a5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3\" # this would allow you to support PDB input files with custom formats, current value is default\n# The remaining options help with the parsing of supported residue information from PDB.\nFMCSC_PDB_R_CONV ???? # you need to change from the default of 1 only for specific issues: 3 is good when you work with files directly from CHARMM, 4 with files directly from AMBER\n# FMCSC_PDB_TOLERANCE_A ???? # increased tolerances (bond angle in degrees) are sometimes to prevent atoms from being rebuilt\n# FMCSC_PDB_TOLERANCE_B ???? ???? # increased tolerances (bond length ratio) are sometimes needed to avoid atoms from being rebuilt") if (pmo.cocho != 2) puts "Are the molecule(s) in the PDB file 'broken' such that all atoms always reside in the central unit cell (GROMACS-style) (y) or are they 'whole' but (possibly) stick out of the unit cell (n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_XYZ_FORCEBOX 1 # 1 assumes molecules are broken into pieces from multiple images on input") else keyf.puts("FMCSC_XYZ_FORCEBOX 0 # 0 assumes molecules are intact on input (not broken into pieces from multiple images)") end end end else puts "For starting a simulation, CAMPARI can create a wide range initial conformations relying on both structure input (PDB) and biased randomization (excluded volume). What scenario best applies to your case?\n(1): The entire system (all atoms) should be read from PDB, and all present residues are CAMPARI-supported.\n(2): The entire system should be randomized, and all present residues are CAMPARI-supported.\n(3): Whatever possible should be read from an incomplete PDB, and all present residues are CAMPARI-supported.\n(4): The entire system should be randomized, which includes residues not supported natively by CAMPARI.\n(5): The entire system including unsupported residues should be read from a single, complete PDB.\n(6): There is an incomplete PDB for structure reading and a separate PDB for specifying unsupported residue topologies." choice = curated_int_answer(gets.chomp,1,6) keyf.puts("#\n# Settings for how to generate the conformation the system initially resides in\n# Note that not all missing atoms in input PDB files involve relevant degrees of freedom.") case choice when 1 keyf.puts("FMCSC_PDBFILE ???? # if this input file is missing or cannot be read, CAMPARI will default to randomization!\nFMCSC_PDB_READMODE 2 # 2 instructs CAMPARI to assume complete input and use Cartesian coordinates\nFMCSC_PDB_HMODE 1 # no need to rebuild") when 3 keyf.puts("FMCSC_PDBFILE ???? # if this input file is missing or cannot be read, CAMPARI will default to randomization!\nFMCSC_PDB_READMODE 3 # 3 is default and instructs CAMPARI to use whatever possible from a PDB (also see tool 'convert_SEQRES_toseq.sh')") when 4 keyf.puts("FMCSC_PDB_TEMPLATE ???? # CAMPARI expects to find sufficient information for all unsupported residues here (see docu)\n# FMCSC_PDB_READMODE 3 # this keyword is ignored by the template reader!\nFMCSC_PDB_HMODE 1 # if problems with hydrogen positions after read-in are found, try 2 instead") when 5 keyf.puts("FMCSC_PDBFILE ???? # this file will be required when using unsupported residues (co-works as the template)\nFMCSC_PDB_READMODE 2 # 2 instructs CAMPARI to assume complete input and use Cartesian coordinates\nFMCSC_PDB_HMODE 1 # no need to rebuild") when 6 keyf.puts("FMCSC_PDBFILE ???? # if this input file is missing or cannot be read, CAMPARI will default to randomization!\nFMCSC_PDB_TEMPLATE ???? # CAMPARI expects to find sufficient information for unsupported residues here (see docu)\nFMCSC_PDB_READMODE 3 # 3 is default and instructs CAMPARI to use whatever possible from a PDB (also see tool 'convert_SEQRES_toseq.sh')\nFMCSC_PDB_HMODE 1 # if problems with hydrogen positions after read-in are found, try 2 instead") end if (choice != 2) keyf.puts("# FMCSC_PDB_INPUTSTRING \"a6,a5,1x,a4,a1,a3,1x,a1,i4,a1,3x,3f8.3\" # this would allow you to support PDB input files with custom formats, current value is default\n# The remaining options help with the parsing of supported residue information from PDB.\nFMCSC_PDB_R_CONV ???? # you need to change from the default of 1 only for specific issues: 3 is good when you work with files directly from CHARMM, 4 with files directly from AMBER\n# FMCSC_PDB_TOLERANCE_A ???? # increased tolerances (bond angle in degrees) are sometimes to prevent atoms from being rebuilt\n# FMCSC_PDB_TOLERANCE_B ???? ???? # increased tolerances (bond length ratio) are sometimes needed to avoid atoms from being rebuilt") if (pmo.cocho != 2) puts "Are the molecule(s) in the PDB file 'broken' such that all atoms always reside in the central unit cell (GROMACS-style) (y) or are they 'whole' but (possibly) stick out of the unit cell (n)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_XYZ_FORCEBOX 1 # 1 assumes molecules are broken into pieces from multiple images on input") else keyf.puts("FMCSC_XYZ_FORCEBOX 0 # 0 assumes molecules are intact on input (not broken into pieces from multiple images)") end end end if ((choice == 1) || (choice == 5)) keyf.puts("FMCSC_RANDOMIZE 0 # this option eliminates randomization wherever feasible") else keyf.puts("# Note that randomization is not parallelized (can be slow - depends on size) and skips some stiff restraint potentials.\n# Be particularly careful if you prescribe chemical crosslinks in your input sequence and ask CAMPARI to randomize such a structure.") end if (choice == 2) keyf.puts("FMCSC_RANDOMIZE 1 # 1 is the same as 2 if there is no PDB input file: exhaustive randomization\nFMCSC_RANDOMATTS 1000 # it is usually feasible to spend a little more than the default settings here\nFMCSC_RANDOMTHRESH 1.0 # a mean target residue pair energy in kcal/mol; too small a value will simply exhaust the attempts (i.e., is harmless)") elsif ((choice == 3) || (choice == 4) || (choice == 6)) keyf.puts("# The supplemental randomization operates hierarchically on missing tails, side chains, extra small molecules, etc.\n# If you have a clashy input structure, you can take advantage of this by removing additional bits from the PDB.\nFMCSC_RANDOMIZE 1 # 1 will randomize as necessary given PDB input\nFMCSC_RANDOMATTS 1000 # it is usually feasible to spend a little more than the default setting here\nFMCSC_RANDOMTHRESH 1.0 # a mean target residue pair energy in kcal/mol; too small a value will simply exhaust the attempts (i.e., is harmless)\n# FMCSC_RANDOMBURIAL ???? # if you are solvating a partially read structure with a dense fluid (molecules placed randomly), this keyword is important, but the default of 0.4 is almost always a good choice") end end keyf.puts("# FMCSC_WATER3S_GEOM 1 # see Tutorial 6 for these\n# FMCSC_WATER4S_GEOM 1 # see Tutorial 6 for these\n# FMCSC_WATER5S_GEOM 1 # see Tutorial 6 for these\n") puts "Keywords for initial structure setup have been added.\n\n" end def nest1_posta(keyf,pmo,isa) puts "CAMPARI's post-processing facility works as follows. It extracts features from a (simulated or analyzed) trajectory on-the-fly, and then works on these features, for example, to do structural clustering, construct Markov state models, or generate SAPPHIRE plots.\n" puts "What selection of functionalities best applies to your case?\n(0): Nothing needed\n(1): Just structural clustering\n(2): Structural clustering to construct and analyze a Markov model\n(3): SAPPHIRE plot\n(4): Just a dimensionality reduction technique\n(5): Markov model of a PIGS data set" pmo.posttask = curated_int_answer(gets.chomp,0,5) if (pmo.posttask == 0) return end keyf.puts("#\n# Settings for the post-processing facility (clustering, MSMs, SAPPHIREs, PCA, et al\nFMCSC_CCOLLECT ???? # this is the fundamental step interval for collecting features") puts "And how many snapshots do you plan (roughly) to analyze in total (keep in mind that the product of dimensions and snapshots is relevant for both cost and memory footprint (a typical 'large' value is 1.0e8, e.g., 1.0e7 snapshots with 10 dimensions)?" pmo.netcsnps = curated_int_answer(gets.chomp,10,2000000000) nest3_clufeats(keyf,pmo) nest3_cludred(keyf,pmo) case pmo.posttask when 1 nest3_clualgo(keyf,pmo) keyf.puts("FMCSC_CLAGT_MSM 1 # every clustering provides a graph in STRUCT_CLUSTERING.graphml; this is the assumed connectivity lag") when 2 nest3_clualgo(keyf,pmo) nest3_clumsm(keyf,pmo) when 3 nest2_clupix(keyf,pmo) when 5 nest3_clualgo(keyf,pmo) nest3_clumsm(keyf,pmo) end end def nest3_cludred(keyf,pmo) puts "CAMPARI's dimensionaliy reduction algorithms involve: a possible transformation of the data (at least centering); computation of a linear transform; possibly rewriting of features. What do you require? Note that you should pick '0' if you did not link LAPACK.\n(0): Nothing\n(1): Principal component analysis without transformation and dimensionality reduction (analysis)\n(2): Time-lagged independent component analysis (tICA) without transformation and dimensionality reduction (analysis)\n(3): Principal component analysis with transformation and possible dimensionality reduction\n(4): Time-lagged independent component analysis (tICA) with transformation and possible dimensionality reduction" pmo.cdimred = curated_int_answer(gets.chomp,0,4) case pmo.cdimred when 1 keyf.puts("# Dimensionality reduction technique settings\nFMCSC_PCAMODE 2 # 2 asks CAMPARI to perform PCA and only report the PCs (not the data)") when 2 keyf.puts("# Dimensionality reduction technique settings\nFMCSC_PCAMODE 4 # 4 asks CAMPARI to perform tICA and only report the tICs (not the data)") when 3 keyf.puts("# Dimensionality reduction technique settings\nFMCSC_PCAMODE 3 # 3 asks CAMPARI to perform PCA, report the PCs, and dump the transformed data\nFMCSC_CREDUCEDIM ???? # a value of 0 means to keep the original data; a positive value means to use the transformed data with the chosen number of ranked dimensions") when 4 keyf.puts("# Dimensionality reduction technique settings\nFMCSC_PCAMODE 4 # 4 asks CAMPARI to perform tICA, report the tICs, and dump the transformed data\nFMCSC_CREDUCEDIM ???? # a value of 0 means to keep the original data; a positive value means to use the transformed data with the chosen number of ranked dimensions") end end def nest3_clualgo(keyf,pmo) terms_unselected = Array.new(3){ |i| i+1 } if (pmo.netcsnps > 100000) terms_unselected.delete(2) end puts "The following structural clustering algorithms implemented in CAMPARI are meaningfully available based on your prior choices:" terms_unselected.each do |i| printf("(%2d): ",i) case i when 1 puts "Tree-based clustering (a density-based algorithm that is extremely efficient)" when 2 puts "Hierarchical clustering (a very general but expensive technique)" when 3 puts "Leader clustering (a very simple technique)" when 4 puts "File-based clustering (preexisting clustering)" end end choice = curated_int_answer(gets.chomp,1,3) keyf.puts("# Clustering algorithm settings") case pmo.cdistance when 1 if (pmo.hfourier == 1) keyf.puts("FMCSC_CRADIUS ???? # the fundamental size threshold for clustering, common guesses would be ~0.2") else keyf.puts("FMCSC_CRADIUS ???? # the fundamental size threshold in degrees for clustering, common guesses would be ~25") end when 2,3 keyf.puts("FMCSC_CRADIUS ???? # the fundamental size threshold in Angstrom for clustering, common guesses would be ~2A") end case choice when 1 keyf.puts("FMCSC_CMODE 5 # 5 asks CAMPARI to use the tree-based algorithm") keyf.puts("FMCSC_BIRCHHEIGHT 16 # the height of the hierarchical clustering tree (can be varied)\nFMCSC_BIRCHMULTI 15 # recluster in second pass al lrelevant levels (can be varied)") if ((pmo.cdistance == 1) && (pmo.hfourier == 1)) keyf.print("FMCSC_CMAXRAD 0.9 # guess for the unitless expected upper range for the inter-snapshot distance spectrum\n") elsif ((pmo.cdistance == 1) && (pmo.hfourier == 0)) keyf.print("FMCSC_CMAXRAD 100.0 # guess for the expected upper range for the inter-snapshot distance spectrum in degrees\n") else keyf.print("FMCSC_CMAXRAD 10.0 # guess for the expected upper range for the inter-snapshot distance spectrum in Angstrom\n") end keyf.print("# FMCSC_CLEADER 1 # this keyword is only of minor relevance in tree-based clustering\nFMCSC_CRESORT 1 # resorting makes sense for analysis consistency but has no actual impact on the results\nFMCSC_CREFINE 0 # refinement is not recommended in tree-based clustering\n") if ((pmo.pmode == 1) || (pmo.pmode == 3)) keyf.puts("FMCSC_BIRCHCHUNKSZ 1 # 1 disables the OpenMP divide-and-merge strategy in the first phase of the tree-based clustering\n# FMCSC_CMERGEDIAM ???? # if the divide-and-merge is active, this would be the key control parameter") end when 2 keyf.puts("FMCSC_CMODE 3 # 3 selects hierarchical clustering: not usable for large data sets\nFMCSC_CLINKAGE 3 # selects mean linkage, generally the most appropriate choice (1/2 are max/min linkage)\nFMCSC_CRESORT 0 # 1 would enable resorting of same-size clusters by centroid index") case pmo.cdistance when 1 if (pmo.hfourier == 1) keyf.puts("FMCSC_CMAXRAD ???? # a threshold for a preliminary clustering (safe heuristic), common guesses would be ~0.35") else keyf.puts("FMCSC_CMAXRAD ???? # a threshold in degrees for a preliminary clustering (safe heuristic), common guesses would be ~40") end when 2,3 keyf.puts("FMCSC_CMAXRAD ???? # a threshold in Angstrom for a preliminary clustering (safe heuristic), common guesses would be ~5A") end keyf.puts("FMCSC_CCUTOFF ???? # the cutoff for the neighbor list in Angstrom; difficult to guess, usually more similar to FMCSC_CMAXRAD than to FMCSC_CRADIUS\n# FMCSC_NBLFILE ???? # if a prior neighbor list file has been generated, it can be reused") when 3 keyf.puts("FMCSC_CMODE 1 # 1 selects the fixed-centroid Leader algorithm\nFMCSC_CLEADER 1 # search points forward but clusters backward (3 is an alternative)\nFMCSC_CRESORT 0 # 1 would enable resorting of same-size clusters by centroid index\nFMCSC_CREFINE 0 # you could set this to 1 for an often inefficient refinement\n") when 4 keyf.puts("FMCSC_CMODE 6 # 6 is file clustering with data read-in, 7 is without data read-in\nFMCSC_CLUFILE ???? # the input file, which is analogous to STRUCT_CLUSTERING.clu\nFMCSC_CLUFILECOL 1 # 1 picks the first column in the file\nFMCSC_CRESORT 0 # 1 would enable resorting of same-size clusters by centroid index\n") end end def nest3_clumsm(keyf,pmo) puts "Some functionalities rely on a linear algebra library, HSL. Is HSL available?" pmo.havehsl = curated_bin_answer(gets.chomp) keyf.puts("# Fundamental keywords\nFMCSC_TMAT_MD 1 # 1 treats only an assumed forward time axis (3 does both)\nFMCSC_INISYNSNAP 0 # some analyses downstream are operating only on a single component or require a reference snapshot: this allows selecting it\n# Transition network inference settings\nFMCSC_CLAGT_MSM ???? # the lag time is fundamental for Markov models - measured in recorded steps (sliding window)") if (pmo.posttask == 5) keyf.puts("# To properly analyze a PIGS simulation, transitions must be rerouted. This relies on the PIGS trace.\n# In a serial or pure OpenMP mode, this will assume a single trajectory as input that is a concatenation of all replicas in order.\nFMCSC_TRACEFILE ???? # the name of the PIGS trace file corresponding to the original run\nFMCSC_RE_TRAJOUT ???? # this must correspond to the original trajectory output interval\nFMCSC_RE_TRAJSKIP ???? # this must correspond to the number of original steps for which trajectory output was suppressed at the beginning\nFMCSC_RE_TRAJTOTAL ???? # this does not refer to steps: it must the number of stored snapshots in each trajectory file\n# For even more complex cases, other generic modifications may be needed in addition.") end keyf.puts("# FMCSC_TRAJBREAKSFILE ???? # by default, CAMPARI assumes linear time connectivity: this can be used to introduce custom breaks (e.g., trajectory ensembles)\n# FMCSC_TRAJLINKSFILE ???? # similarly, this can be used to introduce custom links (e.g., trajectory reseeding)\nFMCSC_TMATREPORT 0 # set to 1 if you want to print a dedicated transition matrix output file\nFMCSC_BRKLNKREPORT 0 # can be used to dump the assumed connectivity to log-output (1 or 2)") puts "Aside from a lag time, a Markov model is influenced by augmentation to fulfill certain constraints, such as detailed balance. Some common options are:\n(0): Leave as is\n(1): Impose detailed balance naively\n(2): Regularize the count matrix with a pseudocount (Bowman) approach\n(3): Split into disconnected components (and analyze separately)\n(4): Split into disconnected components and impose detailed balance separately (naively)\n(5): Split into disconnected components and regularize the count matrix per component with a pseudocount (Bowman) approach\n" choice3 = curated_int_answer(gets.chomp,0,5) case choice3 when 0 keyf.puts("FMCSC_CADDLINKMODE 0 # 0 means that no transitions are added or removed") when 1 keyf.puts("FMCSC_CADDLINKMODE 4 # 4 means that detailed balance is imposed by symmetrization, this will not connect components fully separated to begin with") when 2 keyf.puts("FMCSC_CADDLINKMODE 7 # 7 means that pseudocounts are used to regularize the matrix, destroying its sparsity (use 8 for an alternative)\nFMCSC_CLINKWEIGHT -1 # -1 is a special trigger to set the pseudocount weight to 1/N_states") when 3 keyf.puts("FMCSC_CADDLINKMODE -1 # -1 means that weakly connected components are fully disconnected; no other changes") when 4 keyf.puts("FMCSC_CADDLINKMODE -4 # -4 means that detailed balance is imposed by symmetrization after separating weakly connected components completely") when 5 keyf.puts("FMCSC_CADDLINKMODE -7 # -7 means that weakly connected components are fully split before pseudocounts are used to regularize the matrices separately, destroying their sparsity (use -8 for an alternative)\nFMCSC_CLINKWEIGHT -1 # -1 is a special trigger to set the pseudocount weight to 1/N_states") end if ((choice3 != 2) && (choice3 != 5)) keyf.puts("# FMCSC_CLINKWEIGHT ???? # comes into play only for some options of FMCSC_CADDLINKMODE, especially 7 and 8\n") end keyf.puts("# Steady state and other analysis settings\nFMCSC_CREWEIGHT 1 # 1 means that the steady state of the Markov model is computed and used downstream\n# FMCSC_MAXTIME_ITERS ???? # you could limit the runtime of all MSM iterative algorithms (see docu, in s)") if (pmo.havehsl == 1) puts "Spectral analysis allows, for example, the computation of (multiple) implied time scales. It is based on an eigen decomposition of the transition matrix. What description fits your use case best?\n(0): Nothing needed\n(1): Just analyze implied time scales\n(2): Just calculate steady state populations\n(3): Something else" choice = curated_int_answer(gets.chomp,0,3) if (choice > 0) keyf.puts("# Spectral analysis settings") end case choice when 1 keyf.puts("FMCSC_EIGVAL_MD 2 # 2 means to identify eigenvalues with largest real components\nFMCSC_NEIGV ???? # the second eigenvalue is the slowest time scale; this sets the total number (ordered)\nFMCSC_DOEIGVECT 0 # disables computation and print-out of the corresponding vectors") when 2 keyf.puts("FMCSC_EIGVAL_MD 2 # 2 means to identify eigenvalues with largest real components\nFMCSC_NEIGV 1 # the steady state is the first eigenvector\nFMCSC_DOEIGVECT 1 # ensure computation and print-out of the corresponding vector") when 3 keyf.puts("FMCSC_EIGVAL_MD ???? # 1/2/3 means to identify eigenvalues with largest absolute/real/complex components\nFMCSC_NEIGV ???? # the steady state is the first eigenvector, slower ones are interpreted as time scale information\nFMCSC_DOEIGVECT 1 # ensure computation and print-out of the corresponding eigenvectors") end if (choice > 0) keyf.puts("# FMCSC_EIGTOL ???? # loosening the tolerance can theoretically help with numerically difficult cases in the HSL Arnoldi solver\n# FMCSC_NEIGBLOCKS ???? # can be used to customize the HSL Arnoldi solver (w/ FMCSC_NEIGSTEPS)\n# FMCSC_NEIGSTEPS ???? # can be used to customize the HSL Arnoldi solver (w/ FMCSC_NEIGBLOCKS)\n# FMCSC_NEIGRST ???? # a heuristic restart attempt number for the HSL Arnoldi solver") end puts "Committor probabilities are used to analyze the probability flux in a Markov model, including decompositions into transition paths. Is this something you require (y/n)?" choice2 = curated_bin_answer(gets.chomp) if (choice2 == 0) keyf.puts("FMCSC_DOPFOLD 0 # 0 means to not compute any committor probabilities") else keyf.puts("# Committor probability analyses\nFMCSC_DOPFOLD 1 # 1 tells CAMPARI to compute forward committor probabilities\nFMCSC_CLUFOLDFILE ???? # file to define the target set of states (snapshot-based)\nFMCSC_CLUUNFOLDFILE ???? # file to defined the source (alternative) set of states (snapshot-based)\nFMCSC_PFOLDREPORT 0 # 1 would produce debugging output") if ((choice3 == 0) || (choice3 == 3)) keyf.puts("FMCSC_DOPFOLD_MINUS 1 # if detailed balance is not imposed by FMCSC_CADDLINKMODE, the backwards committor must be calculated explicitly\nFMCSC_CMSMCFEP 10 # 10 asks CAMPARI to use committor probabilties to construct cut-based pseudo free energy profiles (1 would be MFPT)") else keyf.puts("FMCSC_DOPFOLD_MINUS 0 # if detailed balance is imposed by FMCSC_CADDLINKMODE, this computation would be redundant\nFMCSC_CMSMCFEP 8 # 8 asks CAMPARI to use forward committor probabilties to construct a cut-based pseudo free energy profile (1 would be MFPT)\nFMCSC_TPTMODE 2 # classic literature approach (all selectable approaches have drawbacks)\nFMCSC_TPTMAXPATHS 1000 # such a small value is only useful for Markov models with very few states\nFMCSC_TPTMAXFRACTION 0.99 # never reached in practice for real systems\nFMCSC_MAXTIME_TPT 1.0 # a more meaningful termination condition for the algorithm via runtime in h\n# FMCSC_TPTBTNFRACTION ???? # relevant in other choices for TPTMODE\n# FMCSC_TPTMAXKSP ???? # relevant in other choices for TPTMODE") end end else keyf.puts("FMCSC_CMSMCFEP 1 # 1: calculate mean-first passage times to a reference node (see FMCSC_INISYNSNAP) and get cut-based pseudo-FEP") end keyf.puts("# MSMs can be used to generate synthetic trajectories at the state level.\nFMCSC_SYNTRAJ_MD 0 # change this to 1-3 to generate stochastic trajectories at the network state level\n# FMCSC_NSYNTRAJS ???? # this would be the target number of trajectories\n# FMCSC_NSYNSNAPS ???? # this would be the number of steps per trajectory for some modes\n# FMCSC_SYNTRAJOUT ???? # output frequency for these state trajectories\n# FMCSC_ENDSYNSNAP ???? # for some modes, this would define the target end state (snapshot-based)\n# Mean-first passage times are defined between all pairs of states and for small enough numbers of states the full matrix is computable\n# FMCSC_MFPT_MATRIX 1 # requires LAPACK and either CREWEIGHT 1 or presence of HSL\n") end def nest1_parastype(keyf,pmo) case pmo.mpi_mode when 1 keyf.puts("#\n# Parallel simulation settings") puts "\nIn MPI simulation modes with identical simulation conditions for every replica, there are currently two modes of operation:\n(1) Use a simple MPI averaging calculation (multiple independent replicas with data pooling)\n(2) Use Progress Index-Guided Sampling (PIGS), an enhanced sampling technique relying on clever restarts" choice = curated_int_answer(gets.chomp,1,2) case choice when 1 keyf.puts("FMCSC_MPIAVG_XYZ 1 # keeping any possible trajectories separate is generally easier to manage later on") when 2 pmo.pigsflag = 1 puts "How many copies (replicas) do you intend to use (at least 2)?" choice2 = curated_int_answer(gets.chomp,2,100001) gpv = (choice2/2).to_i keyf.print("FMCSC_MPI_PIGS 1 # this enables the PIGS method if FMCSC_MPIAVG is 1\nFMCSC_MPI_GOODPIGS ",gpv," # the number of top-ranked copies to form the reseeding pool (can be changed)\n") puts "And how many steps should there be between reseeding intervals (the total number of reseeding cycles is typically in the 100s or higher)?" choice3 = curated_int_answer(gets.chomp,5,2000000000) keyf.print("FMCSC_REFREQ ",choice3," # the deterministic interval for performing reseeding cycles\n") cc = [1,(choice2*choice3/100000).to_i].max if (cc <= 5) while ((choice3/cc) > 5) cc = cc + 1 end end if (choice3/cc < 5) cc = (choice3/5).to_i end keyf.print("FMCSC_CCOLLECT ",cc," # the collection interval for the features underlying the PIGS heuristic\n") pmo.netcsnps = (choice2*choice3/cc).to_i nest3_clufeats(keyf.pmo) nest2_clupix(keyf,pmo) end when 2 keyf.puts("#\n# Replica exchange simulation settings (a REX move is essentially a Monte Carlo move in replica space)") puts "A replica exchange simulation uses exchanges between elements in a list of conditions to accelerate sampling. This list is specified by an input file you have to write.\n How do you want exchanges of conditions to occur?\n(1): Across all replicas (watch out for diminishing returns)\n(2): Only between neighboring replicas (the standard protocol)?" choice = curated_int_answer(gets.chomp,1,2) keyf.puts("FMCSC_REFILE ???? # every REX-MPI run requires a custom input file to indicate conditions per replica\nFMCSC_REDIM ???? # this has to be the actual number of dimensions of the conditions (often 1) - has to match the input file\n# FMCSC_RETRACE 1 # there is usually no reason to disable the writing of any possible trace file\n") case choice when 1 keyf.puts("FMCSC_RENBMODE 1 # all-to-all exchanges allowed \nFMCSC_RESWAPS ???? # the number of swap attempts per cycle: the number of unique, eligible pairs is N*(N-1)/2\nFMCSC_REFREQ ???? # the step interval for attempting swap cycles (deterministic frequency), values larger than FMCSC_NRSTEPS disable exchanges") when 2 keyf.puts("FMCSC_RENBMODE 2 # only neighbor exchanges allowed: this is cheaper and usually robust\nFMCSC_RESWAPS ???? # the number of swap attempts per cycle: the number of unique, eligible pairs is (N-1)\nFMCSC_REFREQ ???? # the step interval for attempting swap cycles (deterministic frequency), values larger than FMCSC_NRSTEPS disable exchanges") end if (pmo.simtype != 1) keyf.puts("FMCSC_RE_VELMODE 2 # keep directions but scale velocities for target condition as needed (questionable for Hamiltonian replica exchange)") else keyf.puts("FMCSC_REMC_DOXYZ 1 # exchange Cartesian rather than internal coordinates; if communication is an issue, internals should be faster (set to 0)") end end if (pmo.mpi_mode > 0) puts "Keywords specifying the MPI parallel mulitcopy simulation mode have been added." end end def nest3_clufeats(keyf,pmo) if (pmo.pigsflag == 1) puts "In PIGS, the approximate progress index method is used to inform reseeding decisions. This requires defining a metric, the essential hyperparameter of the method." # else # puts "The structural clustering facility is primarily a post-processing operation. During the run, (preliminary) features are extracted from the trajectory." end puts "What kind of coordinates do you want to base your metric on?\n(1): Dihedral angles\n(2): Interatomic distances\n(3): Absolute coordinates (RMSD)" pmo.cdistance = curated_int_answer(gets.chomp,1,3) pmo.hfourier = 0 samealign = 1 pmo.prepmode = 0 wmode = 0 puts "Do you want to use data-derived weights for individual features (calculated by CAMPARI) to shape the metric further (y/n)?" hweights = curated_bin_answer(gets.chomp) keyf.puts("# Feature selection and processing settings\nFMCSC_CFILE ???? # this is an essential input file for clustering to select the degrees of freedom: a complete set will be used if absent") if (pmo.cdistance == 3) if (hweights == 0) puts "Do you want to use the same set for alignment and for computing the metric in the RMSD calculation (y/n)?" samealign = curated_bin_answer(gets.chomp) if (samealign == 1) keyf.puts("FMCSC_CDISTANCE 5 # 5 selects absolute coordinates as source features (RMSD)\nFMCSC_CALIGN 1 # this precedes every distance evaluation with pairwise alignment (not cheap)\nFMCSC_CPREPMODE 0 # none; the only supported option would be 1 (data centering)") else keyf.puts("FMCSC_CDISTANCE 6 # 6 selects absolute coordinates as source features with split sets (RMSD)\nFMCSC_CALIGN 1 # FMCSC_CDISTACE 6 implies that every distance evaluation performs pairwise alignment (not cheap)\nFMCSC_ALIGNFILE ???? # the required input file to define the alignment coordinates (separate from FMCSC_CFILE)\nFMCSC_CPREPMODE 0 # the only supported option would be 1 (data centering)") end else keyf.puts("FMCSC_CDISTANCE 10 # 10 selects absolute coordinates as source features with weights (RMSD)\nFMCSC_CALIGN 0 # weighted RMSD implies that pairwise alignment is not possible (use prealigned data instead)") end elsif (pmo.cdistance == 2) if (hweights == 0) keyf.puts("FMCSC_CDISTANCE 7 # 7 selects interatomic distances as source features") else keyf.puts("FMCSC_CDISTANCE 9 # 9 selects interatomic distances as source features with weights") end keyf.puts("FMCSC_CDISTRANSFORM 0 # you could change this to transform distances, e.g., in a contact map or DRID-like (reciprocal) sense\n# FMCSC_CDISTRANS_P1 ???? # required parameter for distance transforms\n# FMCSC_CDISTRANS_P2 ???? # additional parameter for some distance transforms") nest4_cpreproc(keyf,pmo) elsif (pmo.cdistance == 1) puts "Are the dihedral angles to be represented by their Fourier terms (y/n)?" pmo.hfourier = curated_bin_answer(gets.chomp) if ((hweights == 0) && (pmo.hfourier == 0)) keyf.puts("FMCSC_CDISTANCE 1 # 1 selects raw dihedral angles: these are periodic quantities, which limits some downstream functionalities and introduces ambiguities") elsif ((hweights == 1) && (pmo.hfourier == 0)) keyf.puts("FMCSC_CDISTANCE 2 # 2 selects weighted dihedral angles: these are periodic quantities, which limits some downstream functionalities and introduces ambiguities") elsif ((hweights == 0) && (pmo.hfourier == 1)) keyf.puts("FMCSC_CDISTANCE 3 # 3 selects dihedral angles as sine/cosine terms; this doubles dimensionality and implies circularly correlated variables") elsif ((hweights == 1) && (pmo.hfourier == 1)) keyf.puts("FMCSC_CDISTANCE 4 # 4 selects weighted sine/cosine terms; this doubles dimensionality and implies circularly correlated variables (weights are per dihedral angle)") end if (pmo.hfourier == 0) keyf.puts("FMCSC_CPREPMODE 0 # explicit dihedral angles do not support preprocessing") else keyf.puts("FMCSC_CPREPMODE 0 # sine/cosine terms have finite range and fixed range anyway, which means that options 0 and 3 are most applicable") end end if (hweights == 1) puts "A number of ways of weighting features are implemented. The most common are:\n(0): Leave the original (most often meaningless)\n(1): Weight by inverse, windowed, mean square fluctuation\n(2): Weight by autocorrelation at fixed lag time\n(4): Weight by crossings of the global feature mean\n(8): Weight by crossing of peaks in the feature distributions\n(3), (5), (6), (7), (9) are combined measures (see documentation)" wmode = curated_int_answer(gets.chomp,0,9) if (wmode == 0) keyf.print("FMCSC_CMODWEIGHTS 0 # default weights are not altered: the only nonuniform default weights are those for dihedral angles and based on conformation\n") else keyf.print("FMCSC_CMODWEIGHTS ",wmode," # this selects the type of weights to use: see documentation\n") end if ((wmode == 2) || (wmode == 6) || (wmode == 7)) keyf.print("FMCSC_CLAGTIME ???? # this sets, in number of analyzed snapshots, the fixed lag for autocorrelation-based weights\n") end if ((wmode != 2) || ((wmode == 0) && (pmo.cdistance == 1))) keyf.puts("FMCSC_CWCOMBINATION 0 # how to combine locally adaptive weights per feature (0 is the geometric mean)\n") end if ((wmode != 2) && (wmode != 0)) keyf.puts("FMCSC_CWINDOWSIZE ???? # an essential parameter for most locally adaptive weights: large values make the weights increasingly static") end if (wmode >= 4) keyf.puts("FMCSC_CTRANSBUF ???? # weights based on crossing counts can use an offset to buffer very small values") end if ((wmode == 9) || (wmode == 5) || (wmode == 7)) if (pmo.prepmode >= 3) keyf.puts("# FMCSC_CSMOOTHORDER ???? # this keyword controls the smoothing interval for weights generation but data were already smoothed (FMCSC_CPREPMODE): see above!") else keyf.puts("FMCSC_CSMOOTHORDER ???? # this keyword controls the smoothing size (B-splines) for features purely for the purpose of deriving weights") end end else keyf.puts("# FMCSC_CMODWEIGHTS 0 # if there are no weights in use, this keyword is redundant") end if (pmo.pigsflag == 0) keyf.puts("# FMCSC_CDUMP 0 # turning this to 1 instructs CAMPARI to dump the collected and processed features to file") end end def nest2_clupix(keyf,pmo) keyf.puts("#\n# Clustering settings (note that suggested thresholds assume that no variance normalization or dimensionality reduction of features has occured)\nFMCSC_CMODE 4 # 4 selects progress index mode\nFMCSC_CPROGINDMODE 2 # using the approximate version implies a preliminary tree-based clustering\nFMCSC_BIRCHHEIGHT 12 # the height of the hierarchical clustering tree (can be varied)\nFMCSC_BIRCHMULTI 6 # recluster in second pass the finest 6 levels (can be varied)") if ((pmo.cdistance == 1) && (pmo.hfourier == 1)) keyf.print("FMCSC_CRADIUS ???? # the unitless fine threshold for sin/cos of selected dihedral angles; a common first guess would be 0.2-0.3\nFMCSC_CMAXRAD 0.9 # the unitless expected upper range for the inter-snapshot distance spectrum\n") elsif ((pmo.cdistance == 1) && (pmo.hfourier == 0)) keyf.print("FMCSC_CRADIUS ???? # the unitless fine threshold in degrees for selected dihedral angles; a common first guess would be 25-35\nFMCSC_CMAXRAD 100.0 # the expected upper range for the inter-snapshot distance spectrum in degrees\n") else keyf.print("FMCSC_CRADIUS ???? # the unitless fine threshold in Angstrom; a common first guess would be 1.5-2.5\nFMCSC_CMAXRAD 10.0 # the expected upper range for the inter-snapshot distance spectrum in Angstrom\n") end keyf.print("# FMCSC_CLEADER 1 # this keyword is only of minor relevance in tree-based clustering\nFMCSC_CRESORT ",1-pmo.pigsflag," # resorting makes sense for analysis consistency but has no actual impact on the results\nFMCSC_CREFINE 0 # refinement is not recommended in tree-based clustering\n") if ((pmo.pmode == 1) || (pmo.pmode == 3)) keyf.puts("FMCSC_BIRCHCHUNKSZ 1 # 1 disables the OpenMP divide-and-merge strategy in the first phase of the tree-based clustering\n# FMCSC_CMERGEDIAM ???? # if the divide-and-merge is active, this would be the key control parameter") end keyf.print("#\n# Progress index settings\nFMCSC_CPROGINDSTART -1 # this uses the centroid representative of the largest cluster as the starting point: can be customized\n") if (pmo.pigsflag == 1) keyf.print("FMCSC_CPROGINDRMAX 100 # the key parameter of the approximation: the number of search attempts per point and stage (should be constant and matched roughly to leaf cluster sizes)\nFMCSC_CPROGRDEPTH 3 # search up to 3 extra levels toward root from current starting level; smaller is faster and more approximate\nFMCSC_CPROGMSTFOLD 0 # the technique controlled by this keyword is likely counterproductive for PIGS\n") else keyf.print("FMCSC_CPROGINDRMAX 500 # the key parameter of the approximation: the number of search attempts per point and stage (should be constant and matched roughly to leaf cluster sizes)\nFMCSC_CPROGINDWIDTH ",[100,(pmo.netcsnps/20).to_i].max," # an output-relevant value that should be matched to the mean basin size to be resolved in number of snapshots\nFMCSC_CPROGMSTFOLD 3 # a technique to achieve better kinetic homogeneity in the PIX by outlier reassignment (see docu)\nFMCSC_CPROGRDEPTH 11 # search up to 11 extra levels toward root from current starting level; smaller is faster but more approximate\n# FMCSC_CBASINMAX ???? # would come into play should FMCSC_CPROGINDSTART be 0\n# FMCSC_CBASINMIN ???? # would come into play should FMCSC_CPROGINDSTART be 0\n# FMCSC_CLAGT_MSM ???? # the cut function can use a lag time (default is 1 snapshot); this keyword is more central to downstream analyses\n") end if ((pmo.pmode == 1) || (pmo.pmode == 3)) keyf.print("FMCSC_CPROGRDBTSZ 50 # a setting that helps with threads-parallel performance\n") end end def nest4_cpreproc(keyf,pmo) puts "The following options are available for preprocessing:\n(0): None\n(1): Data are centered\n(2): Data are centered and normalized to unit variance (Z-scored)\n(3): Data are smoothed (using B-splines)\n(4): Data are centered and smoothed (using B-splines)\n(5): Data are centered and normalized to unit variance (Z-scored) and smoothed (using B-splines)" pmo.prepmode = curated_int_answer(gets.chomp,0,5) case pmo.prepmode when 0 keyf.puts("FMCSC_CPREPMODE 0 # 0 means no preprocessing") when 1 keyf.puts("FMCSC_CPREPMODE 1 # 1 means centering (all features have zero means)") when 2 keyf.puts("FMCSC_CPREPMODE 2 # 2 means centering and variance normalization (features become Z-scores)") when 3 keyf.puts("FMCSC_CPREPMODE 3 # 3 means just smoothing (this is a significant alteration of the data)\nFMCSC_CSMOOTHORDER ???? # the number of snapshots corresponding to the B-spline smoothing window; use FMCSC_CDUMP to look at the smoothed features directly") when 4 keyf.puts("FMCSC_CPREPMODE 4 # 4 means centering and smoothing (this is a significant alteration of the data)\nFMCSC_CSMOOTHORDER ???? # the number of snapshots corresponding to the B-spline smoothing window; use FMCSC_CDUMP to look at the processed features directly") when 5 keyf.puts("FMCSC_CPREPMODE 5 # 5 means Z-scoreing and smoothing (this is a significant alteration of the data)\nFMCSC_CSMOOTHORDER ???? # the number of snapshots corresponding to the B-spline smoothing window; use FMCSC_CDUMP to look at the processed features directly") end end def nest2_lrel(keyf,pmo) puts "Since you selected Coulomb interactions, how should their long-range effects be accounted for?" terms_unselected = Array.new(5){ |i| i+1 } if ( (pmo.cocho != 1) || (pmo.stdnbs != 1) || (pmo.mcflag == 1) ) terms_unselected.delete(4) terms_unselected.delete(5) elsif (pmo.ghosting == 1) terms_unselected.delete(5) end choice = 0 terms_unselected.each do |i| printf("(%2d): ",i) case i when 1 puts "Truncation (strict cutoff at specified value)" when 2 puts "Truncation for interactions involving at least one net-neutral group, explicit monopole-based calculation otherwise" when 3 puts "Truncation only for interactions involving two net-neutral groups, explicit atomistic calculation otherwise" when 4 puts "Reaction-field approach wherever possible (implies homogeneous background dielectric)" when 5 puts "Ewald sums wherever possible (implies homogeneous background dielectric)" end end keyf.puts("# Long-range model for Coulomb interactions") choice = curated_int_answer(gets.chomp,1,5) case choice when 1 keyf.puts("FMCSC_LREL_MD 1 # selects truncation for (any possible) dynamics part\n") keyf.puts("FMCSC_LREL_MC 4 # selects truncation for (any possible) Monte Carlo part\n") when 2 keyf.puts("FMCSC_LREL_MD 4 # selects truncation except monopole-monopole for (any possible) dynamics part\n") keyf.puts("FMCSC_LREL_MC 3 # selects truncation except monopole-monopole for (any possible) Monte Carlo part\n") when 3 keyf.puts("FMCSC_LREL_MD 5 # selects truncation only for dipole-dipole for (any possible) dynamics part\n") keyf.puts("FMCSC_LREL_MC 1 # selects truncation only for dipole-dipole for (any possible) Monte Carlo part\n") when 4 keyf.puts("FMCSC_LREL_MD 3 # selects reaction-field method for (any possible) dynamics part\n") keyf.puts("FMCSC_RFMODE 1 # the generalized reaction-field method recognizes the explicit presence of charged species\n") keyf.puts("FMCSC_LREL_MC 3 # selects truncation except monopole-monopole for (any possible) Monte Carlo part (inconsistent with FMCSC_LREL_MD)\n") when 5 keyf.puts("FMCSC_LREL_MD 2 # selects Ewald method for (any possible) dynamics part\n") keyf.puts("FMCSC_EWALD 1 # the particle-mesh Ewald (PME) approach is computationally viable in general; the simple direct Ewald sum approach is not\n") keyf.puts("# FMCSC_EWERFTOL ???? # you could control the accuracy of the tabulation for erfc() assuming it is in use (this is a compile time choice, see installation instructions)\nFMCSC_BSPLINE 6 # the B-spline order for PME\nFMCSC_EWFSPAC 2.0 # sets the Fourier spacing (2.0 corresponds to ~1.0A^3 cubic voxels for charge splining)\n# FMCSC_EWPRM ???? # you could use this parameter to change the relative importance of real- and reciprocal space parts of the infinite sums (you should know what you are doing)\n# FMCSC_EWFFTPLANNER ???? # you can control the plan optimization work conducted by the FFTW library\n# FMCSC_EWSAVEWISDOM ???? # you could use this to save a plan created by expensive optimization\n# FMCSC_EWWISDOMFILE ???? # you could provide an input file with a plan created previously\n") keyf.puts("FMCSC_LREL_MC 3 # selects truncation except monopole-monopole for (any possible) Monte Carlo part (inconsistent with FMCSC_LREL_MD)\n") end if ((choice <= 3) && (choice >= 2)) if (pmo.patching == 1) keyf.puts("# FMCSC_NCPATCHFILE ???? # this specialized patch functionality allows the charge group partitioning / flagging to be modified (see also FMCSC_POLTOL)") end end end def nest2_dofs(keyf,pmo,doholo) pmo.holomode = 0 puts "\nWhat should the degrees of freedom be?\n(1): Internal (torsional and rigid-body) degrees of freedom\n(2): Standard Cartesian coordinates" choice = curated_int_answer(gets.chomp,1,2) case choice when 1 keyf.puts("FMCSC_CARTINT 1 # 1 selects internal coordinates\nFMCSC_TMD_UNKMODE 1 # change this to 3 if you also want to sample otherwise frozen degrees of freedom in supported residues (like CH3 spins)") when 2 keyf.puts("FMCSC_CARTINT 2 # 2 selects Cartesian coordinates (dynamics only)") if (doholo == 1) nest3_holos(keyf,pmo) end end pmo.dofflag = choice end def nest2_relax(keyf,pmo) if (pmo.dofflag == 1) if (pmo.dockrun == 1) puts "\nCAMPARI can try to relax every screened molecule's starting position by a targeted Monte Carlo procedure different from (but possibly similar to) the normal sampler. Do you want to use this (makes most sense in tethered docking with a flexible receptor)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_TMD_RELAX 10.0 # the threshold is not easy to interpret; 0.0 disables this procedure\nFMCSC_TMD_NRRELAX 100 # the number of elementary MC moves per identified residue/dihedral angle to relax") if ((pmo.dckstyle <= 2) && ((pmo.dockteth == 2) || (pmo.dockteth == 3))) puts "\nDo you want to use this technique to prune molecules that are unlikely to fit (useful especially for narrow binding sites and when working with general (not targeted) libraries)?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_MOL2PRUNEMODE 5 # the default is 0 and disables this procedure, be wary of this introducing false negatives") else keyf.puts("FMCSC_MOL2PRUNEMODE 0 # off, which is the default") end end else keyf.puts("# FMCSC_TMD_RELAX 10.0 # the default is 0.0 and disables this procedure\n# FMCSC_TMD_NRRELAX 100 # the number of elementary MC moves per identified residue/entity to relax") end else puts "\nDo you want to try to relax degrees of freedom experiencing high forces in the starting structure by a targeted Monte Carlo procedure? This can be time-consuming but still much more efficient in relaxing specific, isolated clashes in large systems." if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_TMD_RELAX 10.0 # the threshold is not easy to interpret; 0.0 disables this procedure\nFMCSC_TMD_NRRELAX 100 # the number of elementary MC moves per identified residue/entity to relax") else keyf.puts("# FMCSC_TMD_RELAX 10.0 # the default is 0.0 and disables this procedure\n# FMCSC_TMD_NRRELAX 100 # the number of elementary MC moves per identified residue/entity to relax") end end end end def nest2_mini(keyf,pmo) if (pmo.doffl == 1) puts "Note that minimization in internal coordinates is not a trivial procedure (see tutorial #8)." elsif (pmo.doffl == 2) puts "Note that you cannot use holonomic constraints with any minimizer except the stochastic one (under certain conditions)." end addem = 1 pmo.holomode = 0 puts "\nWhat minimization algorithm do you want to use?\n(1): Steepest descent (a safe but inefficient choice)\n(2): Conjugate gradient (generally faster than steepest descent but can get 'stuck' on flat landscapes)\n(3): BFGS (a second order method that is generally fastest in locating nearby minima but can also get 'stuck' on flat landscapes)\n(4): Stochastic minimization (annealing like: uses low-temperature, thermostatted MD in the end)" choice = curated_int_answer(gets.chomp,1,4) case choice when 1 keyf.puts("FMCSC_MINI_MODE 1 # 1 selects steepest descent\n") when 2 keyf.puts("FMCSC_MINI_MODE 2 # 2 selects conjugate gradient\n") when 3 keyf.puts("FMCSC_MINI_MODE 3 # 3 selects the Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS) according to Nocedal\nFMCSC_MINI_MEMORY 10 # the number of previous steps the estimate of the Hessian relies on\n# FMCSC_MINI_UPTOL ???? # you could use this to allow uphill escape routes for BFGS") when 4 keyf.puts("FMCSC_MINI_MODE 4 # 1 selects annealing method\nFMCSC_MINI_SC_HEAT ???? # the fraction of steps used for heating\nFMCSC_MINI_SC_TBATH 0.1 # the temperature in the end in K (should be low)\nFMCSC_ENSEMBLE 1 # the method necessitates the NVT ensemble\nFMCSC_TEMP ???? # this is the initial, hot temperature in K: smaller values increase search locality\nFMCSC_TSTAT 4 # the best general purpose thermostat\nFMCSC_TSTAT_TAU 0.2 # this time constant (in ps) will be modified by the algorithm\nFMCSC_TIMESTEP 0.001 # typical 'safe' integration time step in ps: larger values may be feasible") puts "Do you want to use holonomic constraints with this minimizer?" if (curated_bin_answer(gets.chomp) == 1) addem = 0 keyf.puts("FMCSC_MINI_SC_SDSTEPS 0 # the use of holonomic constraints requires disabling initial SD steps") nest3_holos(keyf,pmo) else keyf.puts("FMCSC_MINI_SC_SDSTEPS ???? # to make the method safer for high energy starting conformations, initial steepest descent steps are offered") end end if (addem == 1) keyf.puts("FMCSC_MINI_GRMS ???? # the convergence criterion for the norm of the gradient: difficult to set a priori, so try very low values first\nFMCSC_MINI_STEPSIZE 1.0 # this controls the initial(!) step size scale for all non-MD minimizers\nFMCSC_MINI_XYZ_STEPSIZE 1.0e-2 # translational step size for all minimizers in Angstrom") if (doffl == 1) keyf.puts("FMCSC_MINI_ROT_STEPSIZE 2.0e-3 # rotational step size in degrees\nFMCSC_MINI_INT_STEPSIZE 5.0e-2 # torsional step size in degrees") end end end def nest2_ensemble(keyf,pmo) puts "\nWhat thermodynamic ensemble should the simulation use?" terms_unselected = Array.new(5){ |i| i+1 } terms_unselected.delete(3) terms_unselected.delete(4) if (pmo.simtype != 1) terms_unselected.delete(5) end if (pmo.simtype != 2) terms_unselected.delete(2) end pmo.ensemble = 0 terms_unselected.each do |i| printf("(%2d): ",i) case i when 1 puts "Canonical (NVT) ensemble" when 2 puts "Microcanonical (NVE) ensemble" when 3 puts "Isothermal-isobaric (NPT) ensemble" when 4 puts "Isoenthalpic-isobaric (NPH) ensemble" when 5 puts "Grand canonical (muVT) ensemble" when 6 puts "Semi-grand canonical ensemble" end end while (pmo.ensemble == 0) pmo.ensemble = curated_int_answer(gets.chomp,1,6) while (!terms_unselected.include?(pmo.ensemble)) puts "Please select an available option." pmo.ensemble = curated_int_answer(gets.chomp,1,6) end end case pmo.ensemble when 1 keyf.puts("FMCSC_ENSEMBLE 1 # NVT\n") when 2 keyf.puts("FMCSC_ENSEMBLE 2 # NVE\n") keyf.puts("FMCSC_TIMESTEP 0.001 # in ps; there is no general value here, but 1fs is a common 'safe' value for molecular systems\n") when 5 keyf.puts("FMCSC_ENSEMBLE 5 # muVT\n") keyf.puts("# Grand-canonical mandatory keywords\nFMCSC_PARTICLEFLUCFILE ??? # the input file for which molecule types to have fluctuate in number\nFMCSC_GRANDMODE 2 # it is recommended to use excess (mode 2) rather than absolute (mode 1) chemical potentials\nFMCSC_PARTICLEFLUCFREQ ??? # this move set frequency is required to be nonzero; make sure you also turn on rigid body moves (FMCSC_RIGIDFREQ and so on)\nFMCSC_NUMCALC 10 # no harm in counting often\n# FMCSC_NUMINST ???? # interval relative to NUMCALC to print instantaneous numbers\nFMCSC_GRANDREPORT 1 # reports are helpful when using custom input files") end end def nest2_coupling(keyf,pmo) if ((pmo.simtype == 3) || (pmo.simtype == 7)) keyf.puts("#\n# Bath coupling and integrator settings") if (pmo.dofflag == 2) if (pmo.holomode > 1) keyf.puts("FMCSC_TIMESTEP 0.002 # 0.002ps is a common safe time step in MD simulations with bond constraints") else keyf.puts("FMCSC_TIMESTEP ???? # e.g., 0.001ps is a common safe time step when sampling H-X bonds explicitly") end else keyf.puts("FMCSC_TIMESTEP ???? # e.g., 0.004ps is a common safe time step when using internal degrees of freedom") end keyf.puts("FMCSC_FRICTION ???? # the friction is usually an artificial parameter: common values are 5/ps to 10/ps") return end if (pmo.holomode > 1) keyf.puts("#\n# Bath coupling and integrator settings\nFMCSC_TSTAT 4 # Bussi-Parrinello is the only real option left due to use of holonomic constraints\nFMCSC_TSTAT_TAU 5.0 # larger coupling times produce more 'natural' dynamics\nFMCSC_TIMESTEP 0.002 # 2fs is a common safe time step in MD simulations with bond constraints") return end keyf.puts("#\n# Bath coupling and integrator settings") if (pmo.dofflag == 1) nest3_tmds(keyf,pmo) end puts "Temperature (or generally bath) coupling is a delicate issue in molecular dynamics. The most common problems are: quenching of fluctuations; equipartition defects; systematic errors. What situation describes your case best?\n(1): Many tightly coupled degrees of freedom (e.g., complex fluid, protein in explicit solvent, etc.)\n(2): Mostly weakly coupled degrees freedom (e.g., implicit solvent simulations in internal coordinates)\n(3): A predominantly constrained system (e.g., in a small molecule screen)\n(4): A mix or something else (e.g., Cartesian coordinates of molecules in an implicit solvent)" choice = curated_int_answer(gets.chomp,1,3) case choice when 1 keyf.puts("FMCSC_TSTAT 4 # Bussi-Parrinello is the way to go in tightly coupled systems (equipartition issues unlikely)\nFMCSC_TSTAT_TAU 5.0 # larger coupling times produce more 'natural' dynamics\nFMCSC_TIMESTEP ???? # e.g., 0.001ps is a common safe time step when sampling H-X bonds explicitly") when 2,3 keyf.puts("FMCSC_TSTAT 2 # Andersen avoids equipartition problems but decorrelates velocities: safe to use primarily in weakly coupled or heavily constrained systems\nFMCSC_TSTAT_TAU 5.0 # larger coupling times produce more 'natural' dynamics\nFMCSC_TIMESTEP ???? # e.g., 0.004ps is a common safe time step when using internal degrees of freedom") when 4 if (pmo.dofflag == 2) keyf.puts("FMCSC_TSTAT 4 # 4 selects Bussi-Parrinello: consider Langevin dynamics as an alternative\nFMCSC_TSTAT_TAU ???? # larger coupling times produce more 'natural' dynamics\nFMCSC_TIMESTEP ???? # e.g., 0.001ps is a common safe time step when sampling H-X bonds explicitly") else keyf.puts("FMCSC_TSTAT ???? # Andersen (2) or Bussi-Parrinello (4) are the only real options\nFMCSC_TSTAT_TAU ???? # larger coupling times produce more 'natural' dynamics\nFMCSC_TIMESTEP ???? # e.g., 0.001ps is a common safe time step when sampling H-X bonds explicitly") end end end def nest2_moveset(keyf,pmo) keyf.puts("# Monte Carlo move set choices") if ((pmo.dockrun == 1) && ((pmo.dckstyle <= 2) || (pmo.dckstyle == 4))) case pmo.dckstyle when 1 keyf.print("FMCSC_RIGIDFREQ 1.0 # all (remaining) moves are of the rigid-body type\nFMCSC_CLURBFREQ 0.0 # only one molecule is moving\n") when 2 keyf.print("FMCSC_RIGIDFREQ 0.25 # fraction of rigid-body moves amongst all moves\nFMCSC_CLURBFREQ 0.0 # only one molecule should be moving\nFMCSC_CRFREQ 0.0 # no concerted rotation moves\nFMCSC_NUCFREQ 0.0 # no nucleic acid backbone moves\nFMCSC_OMEGAFREQ 0.0 # no polypeptide omega moves\nFMCSC_PKRFREQ 0.0 # no pucker moves for proline et al.\nFMCSC_CHIFREQ 0.25 # side chain moves, you might have to set this to 0.0 if the run stops\n# FMCSC_CHIRDFREQ ???? # the default is 10%\n# FMCSC_CHISTEPSZ ???? # the default is 10 degrees\nFMCSC_OTHERFREQ 1.0 # single dihedral angle pivot moves\nFMCSC_OTHERUNKFREQ ???? # the fraction to operate on degrees of freedom in the small molecule (and other unsupporteds, if any); you might have to set this to 1.0 if the run stops\nFMCSC_OTHERNATFREQ ???? # the fraction of remaining OTHER moves to operate on native degrees of freedom (usually redundant with CHIFREQ here)\n# FMCSC_OTHERRDFREQ ???? # the default is 20%\n# FMCSC_OTHERSTEPSZ ???? # the default is 20 degrees\nFMCSC_ALIGN 4 # safe in maintaining absolute reference frame as long as no backbone sampling is attempted\n") when 4 keyf.print("FMCSC_RIGIDFREQ 0.01 # fraction of rigid-body moves amongst all moves (just to avoid boundary artifacts)\nFMCSC_CLURBFREQ 0.0 # only one molecule should be moving\nFMCSC_CRFREQ 0.0 # no concerted rotation moves\nFMCSC_NUCFREQ 0.0 # no nucleic acid backbone moves\nFMCSC_OMEGAFREQ 0.0 # no polypeptide omega moves\nFMCSC_PKRFREQ 0.0 # no pucker moves for proline et al.\nFMCSC_CHIFREQ 0.0 # no side chain moves\nFMCSC_OTHERFREQ 1.0 # single dihedral angle pivot moves\nFMCSC_OTHERUNKFREQ 1.0 # the fraction to operate on degrees of freedom in the small molecule\n# FMCSC_OTHERRDFREQ ???? # the default is 20%\n# FMCSC_OTHERSTEPSZ ???? # the default is 20 degrees\nFMCSC_ALIGN 3\n") end else puts "\nMonte Carlo move sets have, by nature, a large number of parameters. I will try to come up with a reasonable move set based on simple questions." puts "This procedure assumes that you are intending to simulate an atomistic system at standard molecular scales (energy/length).\n" puts "Note that at the move set level, entire classes of degrees of freedom can become frozen. To freeze specific elements within a class of degrees of freedom, resort to custom constraints or preferential sampling strategies instead.\n\n" puts "Q: Is your system dense (are there larger polymers (> 20 residues) that can collapse or formal concentrations of small molecules in excess of 10M)?" if (curated_bin_answer(gets.chomp) == 1) rdf = 0.05 rdt = 0.2 rds = 1.0 stt = 5.0 else rdf = 0.3 rdt = 1.0 rds = 5.0 stt = 20.0 end rdf2 = 2*rdf if (pmo.cocho == 1) rdb = 1.0 else rdb = 1.2 end puts "Q: Should torsional degrees of freedom be frozen (select 'y' if there are none)?" prtro = 1 if (curated_bin_answer(gets.chomp) == 1) keyf.print("FMCSC_RIGIDFREQ 1.0 # all (remaining) moves are of the rigid-body type\n") prtro = 2 else puts "Q: Should rigid-body degrees of freedom be frozen?" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_ALIGN 3 # swiveling rule: options 1-3 are available but 4 would imply de facto rigid-body moves\nFMCSC_RIGIDFREQ 0.0 # disable all rigid-body moves") prtro = 0 else keyf.puts("FMCSC_ALIGN 4 # swiveling rule: this option is the most random and thus versatile (but does not keep reference frames)\n# Rigid-body moves operate on the external degrees of freedom of entire molecules\n# Acceptance rates often reflect strong heterogeneities in the sizes of the molecules (consider preferential sampling to address this)\nFMCSC_RIGIDFREQ ???? # this is close to the root of the MC decision tree, so small values (0.01-0.2) are common\n") end end if (prtro == 1) keyf.print("FMCSC_RIGIDRDFREQ ",rdf2.to_s," # the fraction of moves allowing the molecule to randomize rigid-body coordinates completely\nFMCSC_RIGIDRDBUF ",rdb.to_s," # check documentation for what this is for\nFMCSC_ROTSTEPSZ ",rds.to_s," # rotational step size for local moves in degrees\nFMCSC_TRANSSTEPSZ ",rdt.to_s," # translational step size for local moves in Angstrom\n# FMCSC_COUPLERIGID 0 # setting this to 0 can decouple rotation from translation moves\n# FMCSC_ROTFREQ 0.5 # this would be the subset of pure rotation moves\n# FMCSC_CLURBFREQ ???? # these moves displace several molecules at decreased efficiency: consider using this in dilute systems\n# FMCSC_CLURBMAX ???? # the maximal number of molecules in one of these moves\n# FMCSC_CLURBRDFREQ ???? # the fraction of these moves to apply a completely random perturbation to the group of molecules\n# FMCSC_CLURBROTFREQ ???? # frequency of using subclass of rotational perturbations to groups of molecules\n# FMCSC_CLURBSPINFREQ ???? # frequency within subclass of rotational perturbations to groups of molecules to spin molecules around a parallel axes per molecule\n") end if (prtro == 2) keyf.puts("FMCSC_CHIFREQ 0.0 # sidechain moves: the next major branch off the main trunk of the decision tree\nFMCSC_CRFREQ 0.0 # polypeptide concerted rotation moves are the next major branch\nFMCSC_OMEGAFREQ 0.0 # polypeptide omega moves are the next branch\nFMCSC_NUCFREQ 0.0 # nucleic acid moves are the next major branch\nFMCSC_PKRFREQ 0.0 # polypeptide pucker moves are the next branch\nFMCSC_OTHERFREQ 0.0 # generic single torsion pivot moves are the last major branch (fall-through moves are polypeptide phi/psi moves)") return end puts "Q: What species are there to sample in the simulation?\n(1): Only (poly)peptides and possibly small molecules\n(2): Only nucleic acids and possibly small molecules\n(3): Both (poly)peptides and nucleic acids and possibly small molecules\n(4): Only small molecules\n(5): Only natively unsupported species of no recognized polymer type\n(6): A mix of everything (select this to get a versatile move set for manual customization)" syspoly = curated_int_answer(gets.chomp,1,6) if (syspoly == 5) keyf.print("FMCSC_CHIFREQ 0.0 # sidechain moves: the next major branch off the main trunk of the decision tree\nFMCSC_CRFREQ 0.0 # polypeptide concerted rotation moves are the next major branch\nFMCSC_OMEGAFREQ 0.0 # polypeptide omega moves are the next branch\nFMCSC_NUCFREQ 0.0 # nucleic acid moves are the next major branch\nFMCSC_PKRFREQ 0.0 # polypeptide pucker moves are the next branch\nFMCSC_OTHERFREQ 1.0 # generic single torsion pivot moves are the last major branch (fall-through moves are polypeptide phi/psi moves)\nFMCSC_OTHERUNKFREQ 1.0 # 1.0 means to work exclusively on unsupported residues\nFMCSC_OTHERRDFREQ ",rdf.to_s," fraction to randomize torsion completely\nFMCSC_OTHERSTEPSZ ",stt.to_s," # step size in degrees for local moves of this type\n") return end wantcr = 0 if ((syspoly <= 3) && ((pmo.pmode == 2) || (pmo.pmode == 0))) puts "Q: Do you want to include exact concerted rotation moves (relevant for long, dense polymers)? These moves have significant costs and are not necessarily OpenMP-supported." if (curated_bin_answer(gets.chomp) == 1) wantcr = 1 end end widendof = 0 puts "Q: Do you want to include additional degrees of freedom in native CAMPARI residues for torsional sampling (compare FMCSC_TMD_UNKMODE, options 2-3)? This increases consistency but most of these degrees of freedom are uninteresting (e.g., -CH3 spins).\n" if (curated_bin_answer(gets.chomp) == 1) widendof = 1 end wantcistrans = 0 if ((syspoly != 2) && (syspoly != 4) && (syspoly != 6)) puts "Q: Do you want to allow (poly)peptides to be able to explore the cis/trans isomerization of their omega bonds?\n" if (curated_bin_answer(gets.chomp) == 1) wantcistrans =1 end end want5ring = 0 if ((syspoly != 4) && (syspoly != 6)) puts "Q: Do you want to sample the conformational equilibria of 5-membered flexible rings where supported (proline, sugars)?\n" if (curated_bin_answer(gets.chomp) == 1) want5ring =1 end end if (syspoly == 6) want5ring = 1 wantcistrans = 1 wantcr = 1 widendof = 1 end keyf.print("# Sidechain moves are generic single- or multi-torsion pivot moves applicable to many native CAMPARI degrees of freedom.\n# They are always residue-local and thus have low energy complexity and relatively high acceptance rates.\nFMCSC_CHIFREQ ???? # in general, it is important to have some idea of the importance of sidechain degrees of freedom for the system: common values are 0.1 to 0.3\nFMCSC_CHIRDFREQ ",rdf2.to_s," # full randomization moves are able to jump torsional barriers in, e.g., tyrosine residues\nFMCSC_CHISTEPSZ ",(stt*1.6).to_s," # the step size for local moves in degrees\nFMCSC_NRCHI 2 # the number of sidechain dihedrals to include in a move: 1 is most efficient but eliminates all serendipitous correlation\n") if (wantcr == 1) keyf.puts("# Exact concerted rotation moves solve a loop closure problem across 6 dihedral angles. There are two variants for polypeptides and one for nucleic acids (currently).\n# This is reduced to a root finding problem, which is costly and (currently) not parallelized.\n# These moves are becoming more attractive efficiency-wise the larger the polymers become (relative to pivot moves).\nFMCSC_TORCRMODE 1 # this setting avoids biases caused by poor prerotations\n# FMCSC_UJMAXTRIES ???? # only relevant for FMCSC_TORCRMODE being 2\nFMCSC_UJCRBIAS 10.0 # larger values cause narrower prerotation distributions and lower failure-to-close rates\nFMCSC_UJCRWIDTH 10.0 # larger values (degrees) create wider prerotation distributions and higher failure-to-close to rates\nFMCSC_UJCRSTEPSZ 1.0 # the coarse scanning interval for the gradient-based root search on which to base interval refinement (performance-relevant)") keyf.puts("FMCSC_CRFREQ ???? # exact CR moves are not ergodic on their own: they always need pivot moves as well") if ((syspoly == 1) || (syspoly == 3) || (syspoly == 6)) keyf.puts("FMCSC_ANGCRFREQ 0.0 # these moves would extend the closure degrees of freedom to backbone bond angles -> inconsistent with MC sampler elsewhere\nFMCSC_TORCRFREQ 1.0 # setting this to values smaller than 1.0 would also enable inexact torsional CR moves (discouraged for efficiency reasons)\nFMCSC_TORCROFREQ 0.2 # these are the moves including omega bonds in closure, which tend to be less efficient\nFMCSC_TORCRSCOME 0.1 # this scales the prerotation of omega angles to 10% of the value for phi/psi\nFMCSC_TORCRMIN_DJ 3 # minimum number of pre-rotation dihedrals (omega/phi/psi) for 3-residue closure (phi/psi)\nFMCSC_TORCRMAX_DJ 9 # maximum number of pre-rotation dihedrals for this\nFMCSC_TORCRMIN_DO 3 # minimum number of pre-rotation dihedrals (omega/phi/psi) for 2-residue closure (omega/phi/psi)\nFMCSC_TORCRMAX_DO 9 # maximum number of pre-rotation dihedrals for this\n# FMCSC_CRMODE ???? # only relevant for inexact torsional CR moves\n# FMCSC_CRWIDTH ???? # only relevant for inexact torsional CR moves\n# FMCSC_CRBIAS ???? # only relevant for inexact torsional CR moves\n# FMCSC_CRDOF ???? # only relevant for inexact torsional CR moves\n# FMCSC_UJCRMIN ???? # only relevant for angle-based closures in CR moves\n# FMCSC_UJCRMAX ???? # only relevant for angle-based closures in CR moves\n# FMCSC_UJCRINTERVAL ???? # only relevant for angle-based closures in CR moves\n# FMCSC_UJCRSCANG ???? # only relevant for angle-based closures in CR moves") end if ((syspoly == 2) || (syspoly == 3) || (syspoly == 6)) keyf.puts("FMCSC_NUCCRFREQ 0.2 # note that this move type is not in the same branch (under 'FMCSC_NUCFREQ' instead)\nFMCSC_NUCCRMIN 3 # minimum number of pre-rotation dihedrals\nFMCSC_NUCCRMAX 10 # maximum number of prerotation dihedrals (2 residues): reduce if nucleotide is short\n") end else keyf.puts("FMCSC_CRFREQ 0.0 # 0.0 disables all CR moves for polypeptides\nFMCSC_NUCCRFREQ 0.0 # 0.0 disables the exact concerted rotation branch within the FMCSC_NUCFREQ branch") end if ((syspoly == 1) || (syspoly == 3) || (syspoly == 6)) keyf.print("# Omega moves are highly specialized single torsion pivot moves.\nFMCSC_OMEGAFREQ 0.05 # their overall frequency should be low because the impact of omega angles on conformational equilibria is expected to be low\nFMCSC_OMEGASTEPSZ 5.0 # the step size for local moves in degrees\n") if (wantcistrans == 1) keyf.print("FMCSC_OMEGARDFREQ 0.4 # full randomization moves are able to jump the torsional barrier and isomerize cis/trans\n") else keyf.print("FMCSC_OMEGARDFREQ 0.0 # 0.0 disables full randomization moves, which makes spontaneous cis/trans isomerization very unlikely\n") end else keyf.puts("FMCSC_OMEGAFREQ 0.0 # 0.0 disables all omega pivot moves for polypeptides") end if ((syspoly == 2) || (syspoly == 3) || (syspoly == 6)) keyf.print("# Nucleic acid backbone moves are all moves left after PARTICLEFLUC, RIGID, CHI, CR, OMEGA.\n# If a value for FMCSC_NUCCRFREQ is specified above, it must be less than 1.0 for pivot and sugar pucker moves to occur.\nFMCSC_NUCFREQ ???? # larger values (>0.5) would be more common: for pure NA systems, only OTHER moves matter downstream\nFMCSC_NUCRDFREQ ",rdf.to_s," # the fraction of full randomization moves among pivot backbone moves for NAs\nFMCSC_NUCSTEPSZ ",stt.to_s," # the step size in degrees for local NA backbone pivot moves\nFMCSC_NRNUC 2 # the number of backbone angles to include in a move: 1 is most efficient but eliminates all serendipitous correlation\n") else keyf.puts("FMCSC_NUCFREQ 0.0 # 0.0 disables all nucleic acid-specific moves (pivot, sugar pucker, exact concerted rotation)") end if (want5ring == 1) keyf.print("# Pucker moves for five-membered rings can perturb ring bond angles: with all angles and bonds fixed there are at most 2 solutions.\n# FMCSC_PUCKERSTEP_AN ???? # you could change the local move step size for angles\n# FMCSC_PUCKERSTEP_DI ???? # you could change the local move step size for dihedrals\nFMCSC_SUGARFREQ 0.2 # if possible, enables 5-ring pucker moves for nucleic acids (of all remaining NA moves 0.2 means 20/80 sugar/pivot)\nFMCSC_SUGARRDFREQ 0.2 # fraction of simple (and nonergodic) inversion amongst sugar pucker moves\nFMCSC_PKRFREQ ???? # this value should depend strongly on the fraction of proline or similar residues; note that OTHER and standard phi/psi pivot moves are still further up the tree (meaning FMCSC_PKRFREQ 1.0 would disable those moves)\nFMCSC_PKRRDFREQ 0.2 # fraction of simple (and nonergodic) inversion amongst polypeptide pucker moves\n") else keyf.puts("FMCSC_SUGARFREQ 0.0 # disable 5-ring pucker moves on the sugar groups in nucleic acids\nFMCSC_PKRFREQ 0.0 # disable 5-ring pucker moves on polypeptides") end putothers = 1 keyf.print("# Single torsion pivot moves are the only universal dihedral angle move in CAMPARI (currently).\n# Care has to be taken if certain types of torsions are not meant to be sampled by means of move set selections.\n") if ((widendof == 0) && (wantcistrans == 0)) puts "Q: Are their any unsupported residues with torsions to be sampled in the system?\n" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_OTHERFREQ ???? # the only remaining move class higher in the tree are polypeptide phi/psi pivot moves\nFMCSC_OTHERUNKFREQ 1.0 # all OTHER moves should apply to unsupported residues only") else keyf.puts("FMCSC_OTHERFREQ 0.0 # these are the only moves currently able to sample torsions in unsupported residues") putothers = 0 end elsif ((widendof == 1) && (wantcistrans == 0)) keyf.puts("FMCSC_OTHERFREQ ???? # the only remaining move class higher in the tree are polypeptide phi/psi pivot moves\nFMCSC_OTHERUNKFREQ ???? # this should reflect the proportions of degrees of freedom (unsupported vs supported but normally frozen)\nFMCSC_OTHERNATFREQ 0.0 # 0.0 means to not sample native CAMPARI degrees of freedom (the fall-through is to nonnative dofs in supported residues)") elsif ((widendof == 0) && (wantcistrans == 1)) keyf.puts("FMCSC_OTHERFREQ ???? # the only remaining move class higher in the tree are polypeptide phi/psi pivot moves\nFMCSC_OTHERUNKFREQ ???? # this should reflect the proportions of degrees of freedom (unsupported vs supported but normally frozen)\nFMCSC_OTHERNATFREQ 1.0 # 1.0 means that no nonnative CAMPARI degrees of freedom in supported residues will be sampled") else keyf.puts("FMCSC_OTHERFREQ ???? # the only remaining move class higher in the tree are polypeptide phi/psi pivot moves\nFMCSC_OTHERUNKFREQ ???? # this should reflect the proportions of degrees of freedom (unsupported vs supported but normally frozen)\nFMCSC_OTHERNATFREQ ???? # controls remaining fraction to sample native CAMPARI degrees of freedom (the fall-through is to nonnative dofs in supported residues)") end if (putothers == 1) keyf.print("FMCSC_OTHERRDFREQ ",rdf.to_s," # fraction of these moves to fully randomize the given torsion\nFMCSC_OTHERSTEPSZ ",stt.to_s," # the step size for local moves in degrees\n") end if ((syspoly == 1) || (syspoly == 3) || (syspoly == 6)) keyf.print("# The fallthrough move type of the decision tree are joint phi/psi polypeptide backbone pivot moves.\n# Because of this, these moves do not have their own frequency keyword.\n") keyf.print("FMCSC_PIVOTRDFREQ ",(0.6*rdf).to_s," # fraction of these moves to fully randomize the given torsion\nFMCSC_PIVOTSTEPSZ ",(0.6*stt).to_s," # the step size for local moves in degrees\n") end puts "The default picking probabilities for the degrees of freedom controlled by a given move type are generally uniform (see documentation for details). This behavior can be changed, which includes the possibility to set picking probabilities for individual entries to zero.\nQ: Do you want to sample some degrees of freedom in a class preferentially over others in the same class (requires custom input file)?\n" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("# The preferential sampling facility is a MC-specific way to customize the move set at the level of static picking probabilities.\n# This includes the possibility to introduce degree of freedom-specific constraints.\nFMCSC_PSWFILE ???? # the required custom input file for preferential sampling\nFMCSC_PSWREPORT 1 # reports are useful to enable when using custom file-based inputs") end end puts "Move set keywords have been added." end def nest2_hybrids(keyf,pmo) puts "\nIn a hybrid Monte Carlo / dynamics simulation, alternating segments of each sampler follow one another.\nWhat is the fraction of MC moves ( ]0.0;1.0[ )?" notdone = true mcfrac = curated_num_answer(gets.chomp,0.0,1.0) puts "And what is the longest stretch (in number of elementary steps) you want to consider for any segment?" maxlen = curated_int_answer(gets.chomp,1,pmo.nsteps-1) if (mcfrac > 0.5) mclen = maxlen mdlen = (maxlen/mcfrac - maxlen).to_i else mclen = (maxlen/(1-mcfrac) - maxlen).to_i mdlen = maxlen end keyf.print("# The hybrid cycle settings control the alteration between MC and dynamics segments.\nFMCSC_CYCLE_MC_MAX ",mclen," # maximum stretch length for MC segments\nFMCSC_CYCLE_MC_MIN ",mclen," # smaller values are possible\nFMCSC_CYCLE_MC_FIRST ???? # it can help to have a longer first MC cycle to relax a system and its forces\nFMCSC_CYCLE_DYN_MAX ",mdlen," # maximum stretch length for dynamics segments\nFMCSC_CYCLE_DYN_MIN ",mdlen," # lower values would increase stochasticity but mean that velocity relaxation is harder to control\n") end def nest3_holos(keyf,pmo) puts "\nWhat should the constraints be?\n(1): No constraints (except rigid water molecules)\n(2): All bonds involving hydrogen atoms are constrained\n(3): All bonds are constraints\n(4): All bonds and angles are constrained (will usually fail)\n(5): Constraints are user-customized\n" choice = curated_int_answer(gets.chomp,1,5) if (choice != 1) keyf.puts("FMCSC_SHAKEFROM ???? # where to obtain the 'correct' constraint lengths from\nFMCSC_SETTLEH2O 1 # always settle rigid water models\nFMCSC_SHAKEMETHOD 1 # standard SHAKE (rather than LINCS, which is 4)\nFMCSC_SHAKETOL 1.0e-6 # tight relative tolerance is often feasible\nFMCSC_SHAKEMAXITER 1000 # default\n# FMCSC_LINCSITER 2 # iteration number for rotational shortening if you choose LINCS: use at least 2\n# FMCSC_LINCSORDER 8 # starting expansion order for LINCS if you choose LINCS instead of SHAKE: leave this out in general\n") end case choice when 1 keyf.puts("FMCSC_SHAKESET 1 # 1 disable constraints\n") when 2 keyf.puts("FMCSC_SHAKESET 2 # 2 selects only bonds involving hydrogen\n") when 3 keyf.puts("FMCSC_SHAKESET 3 # 3 selects all bonds\n") when 4 keyf.puts("FMCSC_SHAKESET 5 # 5 selects angle constraints via distances (this fails often with any method/settings)\n") when 5 keyf.puts("FMCSC_SHAKESET 6 # 6 selects custom constraints\nFMCSC_SHAKEFILE ????\n") end pmo.holomode = choice end def nest3_tmds(keyf,pmo) keyf.puts("FMCSC_TMD_INTEGRATOR 1 # the more accurate TMD integrator\nFMCSC_TMD_INT2UP 4 # this number of update steps is usually sufficient\nFMCSC_TMDREPORT 0 # change this to 1 to get a full report on the degrees of freedom\nFMCSC_FUDGE_DYN 1 # this will automatically patch masses at the integrator level so that there are no degrees of freedom with net rotating inertial masses below an internal threshold; if you change this to 0, no mass alteration takes place; the time step will then be limited to ca. 3-4fs unless there are no such degrees of freedom (use FMCSC_TMDREPORT to figure this out)") end def nest3_frzs(keyf,pmo) puts "\nGenerally speaking, do you want to exclude selected degrees of freedom from being sampled by any sampler (custom constraints)?\n" if (curated_bin_answer(gets.chomp) == 1) if (pmo.dofflag == 1) keyf.puts("# Custom constraints provide a way to exclude dihedral angles and rigid body coordinates from the list of sampled degrees of freedom.") if ((pmo.simtype == 2) || (pmo.simtype == 3) || (pmo.simtype == 6)) keyf.puts("# Note that you have to use input mode 'A' to access, e.g., individual cardinal direction for molecular translation or individual torsions in a side chain.") elsif ((pmo.simtype == 5) || (pmo.simtype == 7) || (pmo.simtype == 8)) keyf.puts("# Should you choose input mode 'A', the MC portions will not respect these constraints. This is likely undesirable!") end else keyf.puts("# Custom constraints provide a way to exclude specific cardinal directions of atoms from the list of sampled degrees of freedom.") if (pmo.holomode > 1) keyf.puts("# This is incompatible with holonomic constraints as soon as a (partially) frozen atom is also part of a holonomic constraint group.") end end keyf.puts("FMCSC_FRZFILE ???? # the required input file for reqeusting custom constraints\nFMCSC_FRZREPORT 1 # reports are useful to enable when using custom file-based inputs\n") if (pmo.simtype != 1) keyf.puts("# FMCSC_SKIPFRZ 1 # consider enabling this whenever possible") end end end def nest3_wl(keyf,pmo) puts "\nDo you want to use a flat-histogram method in this Monte Carlo run (you will have to set almost all parameters manually).\n" if (curated_bin_answer(gets.chomp) == 1) keyf.puts("FMCSC_MC_ACCEPT 3 # this activates the Wang-Landau flat histogram methodology\nFMCSC_WL_MODE ???? # reaction coordinate histograms are easier to handle than energy ones\n# FMCSC_WL_RC ???? # might have to be specified\n# FMCSC_WL_MOL ???? # might have to be specified\nFMCSC_WL_BINSZ ???? # a fundamentally important keyword for histogram resolution (static!)\nFMCSC_WL_EXTEND ???? # dynamically extending histograms makes convergence even more challenging\nFMCSC_WL_MAX 1.0 # an (initial) upper bound for the expected values: depends entirely on other choices\n# FMCSC_WL_GINITFILE ???? # it can be easier to create a starting input file with a flat distribution within a regular-space histogram and use FMCSC_WL_EXTEND 1 rather than specify the histogram parameters/rules through keywords\nFMCSC_WL_FREEZE ???? # this parameter is difficult to make a recommendation for\nFMCSC_WL_HUFREQ ???? # must be relatively small (probably 10-100) for anything to proceed\nFMCSC_WL_HVMODE 2 # option 1 is more rigorous but can lead to extremely slow convergence\nFMCSC_WL_FLATCHECK ???? # step interval for how frequently to evaluate whether the convergence parameter f should be reduced\nFMCSC_WL_DEBUG 1 # to have more information from a running WL job") else keyf.puts("FMCSC_MC_ACCEPT 1 # Metropolis criterion") end end #s1 = [*"45671".match('[0-9]*')].join #s2 = "45671" #print s1,s2,String(s1 == s2) #v1 = gets.chomp.is_numeric? #v2 = gets.chomp.is_integer? #print String(v1),String(v2) puts "Welcome to the interactive decision tree for populating CAMPARI key-files with the necessary keywords." puts "Please enter the name for the new key-file:" keyn = gets.chomp #if File::exists?(keyn) # File.delete(keyn) #end keyf = File.new(keyn,"w") keyf.puts("# Mandatory keywords\nPARAMETERS ???? # the parameter file contains the database for (largely energy-related) parameters for CAMPARI-supported systems\nFMCSC_SEQFILE ???? # the sequence file defines the system (as a sequence of residues) to simulate or analyze\n#\n# Fundamental keywords\nFMCSC_UNSAFE 0 # setting this to 1 converts some fatal errors to warnings (often needed when working with unsupported residues)\nFMCSC_BASENAME yourcalc # this is the default name used in output files, change it to customize\nFMCSC_FLUSHTIME 1.0 # in minutes: the time interval for flushing the buffers of most output files, meaning forcing them to be up to date\nFMCSC_PDB_AUXINFO 2 # include extra info in PDB files, which is not the default (0); use 4 to strip the header\nFMCSC_SEQREPORT 1 # the most basic report on system composition: almost always useful but not turned on by default\nFMCSC_UAMODEL 0 # this can be changed to switch from all-atom to united atom models, e.g., GROMOS, CHARMM19, for all supported residues") puts 'Please answer questions with y/n or with numbers (as listed).' #mpi_mode = 0 #pmode = 0 pmo = PMS.new nest1_parallelism(keyf,pmo) pmo.cocho = 0 pmo.simtype = 0 pmo.nsteps = 0 pmo.mcflag = 0 pmo.pigsflag = 0 pmo.dockrun = 0 pmo.docksrun = 0 pmo.dckstyle = 0 pmo.dockteth = 0 pmo.dockrefm = 0 pmo.dofflag = 0 pmo.havetmpl = 0 pmo.netcsnps = 1000000 pmo.ensemble = 1 puts "What general type of run do wish to perform?\n(1): A simulation task for a single system\n(2): A replicated simulation or scoring task for virtual screening / computational docking\n(3): A (particle) trajectory-based analysis task (includes file conversion)" if (pmo.pmode <= 1) puts "(4): A pure (abstract) post-processing task (no implied particle representation)\n" choice = curated_int_answer(gets.chomp,1,4) else choice = curated_int_answer(gets.chomp,1,3) end if (choice == 2) nest1_mol2mode(keyf,pmo) end if (choice < 4) nest1_container(keyf,pmo) end if (choice == 3) keyf.puts("#\n# Enable trajectory analysis mode\nFMCSC_PDBANALYZE 1 # an analysis task proceeds as a simulation only that steps are taken by advancing to the next snapshot\nFMCSC_DYNAMICS 2 # preferring this has technical reasons encountered with big systems") nest1_afiletype(keyf,pmo) nest1_hamiltonian(keyf,pmo,2) nest1_insta(keyf,pmo,1) nest1_analysis(keyf,pmo,1) nest1_posta(keyf,pmo,1) elsif (choice < 3) keyf.puts("#\nFMCSC_PDBANALYZE 0 # the default (0) is to perform a simulation task") nest1_strucinit(keyf,pmo) if (pmo.docksrun == 1) nest1_hamiltonian(keyf,pmo,1) else nest1_stype(keyf,pmo) if ((pmo.dockrun == 1) && ((pmo.dckstyle <= 2) || (pmo.dckstyle == 4))) if (pmo.mpi_mode == 1) puts("Note that it is somewhat unusual to run a small molecule screen in data pooling mode. While it is theoretically possible to run progress index-guided sampling (PIGS) in such a screen, this tool does not support creating the options for this particular combination of methodologies at the moment.\n") end else nest1_parastype(keyf,pmo) end nest1_hamiltonian(keyf,pmo,1) if (pmo.dockrun != 1) nest1_insta(keyf,pmo,0) nest1_analysis(keyf,pmo,0) end end else puts(' ') puts(' ') puts(' /---------\ ') puts(' / \ ') puts(' / \ ') puts(' / \ ') puts(' ! XXXX XXXX ! ') puts(' ! XXXX XXXX ! ') puts(' ! XXX XXX ! ') puts(' ! X ! ') puts(' --\ XXX --- ') puts(' ! ! XXX ! ! ') puts(' ! ! ! ! ') puts(' ! I I I I I ! ') puts(' ! I I I I ! ') puts(' \ / ') puts(' -- -- ') puts(' \---/ ') puts(' XXX XXX ') puts(' XXXX XXXX ') puts(' XXXXX XXXXX ') puts(' XXX XXX ') puts(' XXX XXX ') puts(' XXXXX ') puts(' XXX XXX ') puts(' XXX XXX ') puts(' XXX XXX ') puts(' XXXXX XXXXX ') puts(' XXXX XXXX ') puts(' XXX XXX ') puts(' ') puts(' ') puts(" Unfortunately you have reached a dead end. ") puts(" Hopefully future versions of this tool will ") puts(" address this problem. ") puts(' ') puts(" (the ASCII art is a CHARMM reference) ") return end keyf.close keywords_active = [] keywords_comment = [] final_keyfile = [] uchoices = [] File.readlines(keyn).reverse_each { |s| if ((s.split(" ")[0] == "#") && ((s.split(" ")).length > 1)) if ( ((s.split(" ")[1]).split("_",limit=2))[0] == "FMCSC") # puts ((s.split(" ")[1]).split("_",limit=2))[1] keywords_comment.push(s.split(" ")[1]) end final_keyfile.push(s) elsif ((((s.split(" ")[0]).split("_",limit=2))[0] == "FMCSC") || (s.split(" ")[0] == "PARAMETERS")) if (!keywords_active.include?(s.split(" ")[0])) keywords_active.push(s.split(" ")[0]) if (s.split(" ")[1] == "????") puts s puts "Do you want to put in the above missing specification now (y/n)?\n" if (curated_bin_answer(gets.chomp) == 1) puts "Please type the value (no formatting checks occur):\n" subby = gets.chomp uchoices.push(subby) s.sub!("???? ???? ????",subby) s.sub!("???? ????",subby) s.sub!("????",subby) else uchoices.push((s.split(" ",limit=2)[1]).split("#")[0]) puts "OK, you will have to adjust this manually (dysfunctional key-file as it is)." end end final_keyfile.push(s) end else final_keyfile.push(s) end } longest = [10,((uchoices.sort { |x,y| y.length <=> x.length })[0]).length].max + 2 keyf = File.new(keyn,"w") final_keyfile.reverse_each { |s| if ((((s.split(" ")[0]).split("_",limit=2))[0] == "FMCSC") || (s.split(" ")[0] == "PARAMETERS")) keyf.printf("%-25s %-#{longest}s # %s",s.split(" ")[0],(s.split(" ",limit=2)[1]).split("#")[0],(s.split(" ",limit=2)[1]).split("#",limit=2)[1]) else if ((s.split(" ")[0] == "#") && ((s.split(" ")).length > 1)) if ( ((s.split(" ")[1]).split("_",limit=2))[0] == "FMCSC") keyf.printf("# %-23s %-#{longest}s # %s",s.split(" ")[1],(s.split(" ",limit=3)[2]).split("#")[0],(s.split(" ",limit=3)[2]).split("#",limit=2)[1]) else keyf.puts(s) end else keyf.puts(s) end end } keyf.close puts keywords_active begin # do something rescue # handle exception else # do this if no exception was raised ensure # do this whether or not an exception was raised end