Home

ESPResSo User's Guide

1. V r r the force note the factor 1 r The values of r are assumed to be equally distributed between min and max With a fixed distance of Tmax min Npoints 1 the distance values r in the file are ignored and only included for human readability 5 1 2 Lennard Jones interaction Syntax inter typel type lennard jones O Teut Csnifrelauto rof Tap Tmin Required features LENNARD_JONES Description This command defines the traditional 12 6 Lennard Jones interaction between particles of the types type1 and type2 The potential is defined by 5 1 P A Y Comite if fmin Tot lt T lt Tout Toft 0 otherwise The traditional Lennard Jones potential is the work horse potential of particle particle interactions in coarse grained simulations It is a simple model of the van der Waals interaction and is attractive at large distance but strongly repulsive at short distances rof o corresponds to the sum of the radii of the interaction particles at this radius VLJ r 4ecsnife The minimum of the potential is at r rof 20 At this 44 value of r Vis r e 4ecsnift The attractive part starts beyond this value of r Teut determines the radius where the potential is cut off If conifi is not set or it is set to the string auto the shift will be automatically computed such that the potential is continuous at the cutoff radius If rog is not set it is set
2. 188 inter bondid angle_harmonic K o inter bondid angle_cosine K do inter bondid angle_cossquare K do Required features BOND_ANGLE inter bondid dihedral n K p inter coulomb 0 0 inter coulomb inter coulomb parameters inter coulomb lg p3m feu mesh meshz meshy meshz cao alpha Required features ELECTROSTATICS inter coulomb lg p3m tune tunev2 accuracy accuracy r_cut Test mesh mesh cao cao alpha a Required features ELECTROSTATICS inter coulomb epsilon metallic epsilon n_interpol points mesh_off zoff yoff zoff inter coulomb lg dh K Teut Required features ELECTROSTATICS inter coulomb lg mmm2d mazimal_pairwise_error fixed_far_cutoff dielectric m p dielectric contrasts A Aj Required features ELECTROSTATICS inter coulomb lg mmmid switch_radius bessel_cutoff maximal_pairwise_error inter coulomb lg mmmid tune mazimal_pairwise_error Required features ELECTROSTATICS inter coulomb lg memd f_mass mesh epsilon Required features ELECTROSTATICS inter coulomb elc mazimal_pairwise_error gap_size far_cutoff noneutralization dielectric m dielectric contrasts A Aj Required features ELECTROSTATICS iccp3m n_induced_charges convergence convergence_criterion areas areas normals normals sigmas sigmas epsilons epsilons eps_out eps_out relax relaration_parameter max_iterations maz_iterations ext_field ext_field Required features ELE
3. Heiko Schmitz and Florian Muller Plathe Calculation of the lifetime of positronium in polymers via molecular dynamics simulations J Chem Phys 112 2 1040 1045 2000 doi 10 1063 1 480627 URL http link aip org link JCP 112 1040 8 1 19 J Smiatek M P Allen and F Schmidt Tunable slip boundaries for coarse grained simulations of fluid flow Eur Phys J E 26 115 2008 5 9 1 E R Smith Electrostatic energy in ionic crystals Proc R Soc Lond A 375 475 505 1981 T Soddemann B D nweg and K Kremer A generic computer model for am phiphilic systems Eur Phys J E 6 409 2001 5 1 41 T Soddemann B Dinweg and K Kremer Dissipative particle dynamics A useful thermostat for equilibrium and nonequilibrium molecular dynamics simulations Phys Rev E 68 046702 2003 R Strebel Pieces of software for the Coulombic m body problem Dissertation ETH Z rich 1999 URL http e collection ethbib ethz ch show type diss nr 221 52 53 54 55 56 57 S Succi The lattice Boltzmann equation for fluid dynamics and beyond Oxford University Press USA 2001 A P Thompson S J Plimpton and W Mattson General formulation of pressure and stress tensor for arbitrary many body interaction potentials under periodic boundary conditions Journal of Chemical Physics 131 154107 2009 8 1 21 C Tyagi M Stizen M Sega M Barbosa S Kantorovich and C Holm An itera
4. gencode arch compute_20 code sm_20 will compile code only for Fermi cards Default is to compile for compute models 1 1 and 2 0 i e everything with a G90 chip or newer Note that we require at least compute model 1 1 3 2 make Compiling testing and installing ESPResSo The command make is mainly used to compile the ESPResSo source code but it can do a number of other things The generic syntax of the make command is make options target variable valuel When no target is given the target all is used The following targets are available all Compiles the complete ESPResSo source code The variable myconf can be used to specify the name of the configuration header to be used check Runs the testsuite By default all available tests will be run on 1 2 3 4 6 or 8 processors Which tests are run can be controlled by means of the variable tests which processor numbers are to be used can be controlled via the variable processors Note that depending on your MPI installation MPI jobs can only be run in the queueing system so that ESPResSo will not run from the command line In that case you may not be able to run the testsuite or you have to directly submit the testsuite script testsuite test sh to the queueing system Example make check tests madelung tcl processors 1 2 will run the test madlung tcl on one and two processors clean Deletes all files that were created during the compilation 22 mostlyclean Del
5. observable new needs_radial_profile_specs other parameters center lt cx gt lt cy gt lt cx gt maxr mazr minz minz maxz mazz rbins rbins phibins phibins zbins zbins correlation new obs1 id1 obs2 d2 corr_operation operation dt dt tau_max tau mar tau_lin tau_lin compress1 name compress2 name correlation correlation n_corr correlation id autoupdate start stop correlation id update correlation id finalize correlation id write_to_file filename correlation id print correlation id print average1 variancel correlation_time correlation id print average_errorbars blockfile channel write variable varnamel varname2 blockfile channel write variable all blockfile channel write tclvariable varnamel varname2 blockfile channel write tclvariable all blockfile channel write tclvariable reallyall blockfile channel write particles what range all blockfile channel write bonds range blockfile channel write interactions blockfile channel write random blockfile channel write bit_random blockfile channel write seed blockfile channel write bitseed m j 1 pun Ep E N ER N j E w eN 00 1 00 Ko j Ko NO N N N w m N w N GO E N GO 194 blockfile blockfile blockfile blockfile blockfile blockfile blockfile channel channel channel chan
6. ESPResSo ESPResSo User s Guide for version 3 2 0 May 10 2013 Contents Introduction 1 1 Guiding principles o aa a a 1 2 Available simulation methods a a dr a eek ae lA A ee e TA On UMS i oop os hee ts eo he es hs HR ee ee eee ele Pe eee a Die Site oe ade Oh See ees a Cheah eee e ete GU a A se Ree ee As gs E ces El E ee E a E 2 1 Quick installation oo ee 2 2 Running ESPResSo 2 2 ee 2 3 Creating the first simulation seript o e 24 tutorial tel p ps eee ke sortear BGS bebe Denes a ai IA a ad aa iaa 4 Setting up particles 4 1 part Creating single particles o o 4 2 Creating groups of particle o o eee eee 4 3 constraint Setting up constraints o e 44 Virtual Sites via a a aaa ee es 4 5 Grand canonical feabulel e 5 Setting up interactions aed aes ee eee gees se fe dt a se ee ia oa ee Se Nd De Pees 5 3 Bonded interactions E A ee He Se et eee CO deat tes See keer Se eee ae 5 6 Dihedral interactions 5 7 Coulomb interaction 12 12 13 13 19 20 20 22 23 24 26 26 30 36 39 41 5 8 Dipolar interaction 5 9 Special interaction commands b Setting up the system 6 1 setmd Setting global variables o o 6 2 thermostat Setting up the thermostat 6 3 nemd S
7. GHOST_DEBUG Cellsystem output specific to the handling of ghost cells and the ghost cell communication GHOST_FORCE_DEBUG VERLET_DEBUG Debugging of the Verlet list code of the domain decomposition cell system LATTICE_DEBUG Universal lattice structure debugging HALO_DEBUG GRID_DEBUG PARTICLE_DEBUG Output from the particle handling code P3M_DEBUG ESR_DEBUG debugging of P Ms real space part ESK_DEBUG debugging of PMs k space part EWALD_DEBUG FFT_DEBUG Output from the unified FFT code 201 MAGGS_DEBUG RANDOM_DEBUG FORCE_DEBUG Output from the force calculation loops PTENSOR_DEBUG Output from the pressure tensor calculation loops THERMO_DEBUG Output from the thermostats LJ_DEBUG Output from the Lennard Jones code MORSE_DEBUG Output from the Morse code FENE_DEBUG ONEPART_DEBUG Define to a number of a particle to obtain output on the forces calculated for this particle STAT_DEBUG POLY_DEBUG MOLFORCES_DEBUG LB_DEBUG Output from the lattice Boltzmann code VIRTUAL_SITES_DEBUG ASYNC_BARRIER Introduce a barrier after each asynchronous command completion Helps in detection of mismatching communication FORCE_CORE Causes ESPResSo to try to provoke a core dump when exiting unex pectedly MPI_CORE Causes ESPResSo to try this even with MPI errors 202 C Sample scripts In the directory ESPResSo samples you find several scripts that can serve as samples how to use ESPResSo Ij liquid tcl Simp
8. piston global variable polymer Tcl command prepare_vmd_connection Tcl command pressure principal axis of the moment of inertia 101 quick reference of Tcl commands radial distribution function 110 radial distribution function g r 100 radius of gyration of a chain 08 random number generators 129 random seed Rattle Shake algorithm requirements rigid bond salt Tcl command setmd Tcl command shear rate method skin global variable smooth step interaction soft sphere interaction source directory stored configurations 112 stress tensor 105 stretching force structure factor S q subtracted Lennard Jones bond tabulated bond interactions tabulated interaction Tcl global variables Tcl commands analyze 97 blockfile cellsysten 81 change_volume constraint copy_particles correlation 120 counterions crosslink diamond 33 dielectric fsi 154 iccp3m 67 icosaeder imd integrate inter 43 invalidate_system 1b metadynamics nemd observable 115 parallel_tempering 227 part 26 writepsf Tcl command polymer 30 writevcf Tcl command prepare_vmd_connection writevsf Tcl command salt writevtk Tcl command setmd thermostat 76 uwerr 113 velocities writepdb writepdbfoldchains 135 writepdbfoldtopo 135 writepsf 135 writevcf 133 writevsf 132 writevtk Tcl Tk temperature g
9. 115 115 120 127 127 127 132 134 135 136 138 11 3 Checking for features of ESPResSo 2 o 0 Lattice Boltzmann 12 1 Setting up a LB fluid 2 e e 12 2 LB as a thermostat 12 3 Reading and setting single lattice nOdes 12 4 Visualization 12 5 Setting up boundary conditions 202000000 12 6 Choosing between the GPU and CPU implementations 12 7 Electrohydrodynamics 2 2 0 13 1 Membranes 13 2 Parameters 13 3 Geometry 4 External package mbtools 14 1 Introduction 14 2 Installing and getting started o e e e 14 3 The main tcl Script o e ee 14 4 Analysis 14 5 System generation 14 6 Utils 14 7 mmsgl 5 Under the hood 15 1 Internal particle organization ooo ee O Getting involved a ee root ana a oe a Gis A ER Gee Gen ae oe Pai fh bath eee adh ay aaa ne ee A ESPResSo quick reference B 1 General features e a a B2 Interactions e s ute e wet o BOR OS ed SQ eh ee eB ci B 3 Debug Messages ee C Sample scripts D Maxwell Equations Molecular Dynamics MEMD D 1 Equations of motion 148 148 150 150 151 152 153 153 154 155 155 155 158 158 159 160 162 166 173 179 182 182 184 184 185 185 185 186 197 197 199 201 203 204 204 a
10. A tcl list of sphere radii corresponding to the radii for each colloid type in the system If this is non zero the density profile will be calculated in spherical shells about the colloids in the system identified via colloidmoltypes or if colloidmoltypes is not set then the system center of mass is assumed for the colloid vesicle center e nogrid If this is set a grid mesh will not be used to refine the density profile calculation by taking into account vertical differences between mesh points Calculates the number density of each of the beadtypes given in beadtypes as a function of the vertical distance from the bilayer midplane Lipids are also sorted according to their orientation and assigned to upper or lower leaflets accordingly Thus for a system with 3 beadtypes we would obtain 6 columns of output corresponding to 0 lower 1 lower 2 lower 2 upper 1 upper 0 upper where the number refers to the bead type and upper or lower refers to the bilayer leaflet energy verbose output time_vs_energy Obtains the internal energies of the system from the analyze energy command of ESPResSo flipflop verbose output time_vs_flip Makes a call to the analyze get_lipid_orients command of ESPResSo and compares this with a reference set of lipid orients obtained at the start of the simulation with setup_analysis Based on this comparison the number of lipids which have flipped from their original positions is calculated fluctuat
11. Aj41 9 z Bitp Bi4p 1 9 4 For example in the case of velocity autocorrelation function the above mentioned differ ence has a small value and a random sign i e different contributions cancel each other On the other hand in the of the case of mean square displacement the difference is al ways positive resulting in a non negligible systematic error A more general discussion is presented in Ref 42 126 10 Input Output 10 1 No generic checkpointing One of the most asked for feature that seems to be missing in ESPResSo is checkpointing i e a simple way to tell ESPResSo to store and restore the current state of the simulation and to be able to write this state to or read it from a file This would be most useful to be able to restart a simulation from a specific point in time Unfortunately a simple command checkpoint could not work out of two reasons The main reason is that ESPResSo has no way to determine what information constitutes the actual state of the simulation On the one hand ESPResSo scripts sometimes use Tcl variables that contain essential information about a simulation e g the stored values of an observable that was computed in previous time steps counters etc These would have to be contained in a checkpoint However not all Tcl variables are of interest For example Tcl has a number of automatically set variables that contain information about the hostname the machine type etc These variab
12. Fsaf eS 2 Cee ar q qj a D 12 D 6 For which systems to use the algorithm Although it is not very well known by now this algorithm is a promising alternative to the often used Ewald based methods The main advantages and disadvantages shall be named here However it is still best to understand the concept of the algorithm and figure out for yourself if it may be an option The fields are not calculated for an arbitrary charge distribution but updated from the last solution Therefore particles should not move too much between timesteps less than a lattice cube No procedure for error tuning yet You have to adjust the parameters and deter mine the error yourself Only 3D periodic systems are possible for now With the given interpolation scheme the short range part of the potential is highy underestimated when two particles are in the same lattice cube The initialization routine scales with O N and takes a long time for larger and also inhomogenous systems The algorithm is a local update scheme and spatially varying properties can be applied in the future Because of the locality the algorithm itself scales O N and has a big advantage in speed for larger systems Because of the locality it is highly parallelized It is fast 208 The last item is of course dependent on the system properties But if the charges are evenly distributed and the system is not too sparse this algorithm
13. Taiwan Dec 15 19 2009 ce title Jiri Kolafa and John W Perram Cutoff errors in the ewald summation formulae for point charge systems Molecular Simulation 9 5 351 368 1992 5 7 11 5 7 11 5 8 1 H J Limbach and C Holm Single chain properties of polyelectrolytes in poor solvent J Phys Chem B 107 32 8041 8055 2003 8 1 18 D Magatti and F Ferri Fast multi tau real time software correlator for dynamic light scattering Applied Optics 40 24 4011 4021 AUG 20 2001 ISSN 0003 6935 doi 10 1364 A0 40 004011 A C Maggs and V Rosseto Local simulation algorithms for coulombic interactions Phys Rev Lett 88 196402 2002 5 7 5 D D 4 Bernward A Mann The Swelling Behaviour of Polyelectrolyte Networks PhD thesis Johannes Gutenberg University Mainz Germany December 2005 6 2 4 S Marsili G F Signorini R Chelli M Marchi and P Procacci Orac A molec ular dynamics simulation program to explore free energy surfaces in biomolecular systems at the atomistic level J Comp Chem 31 1106 2009 7 7 B Mehlig D W Heermann and B M Forrest Hybrid monte carlo method for condensed matter systems Phys Rev B 45 679 685 1992 P Nikunen M Karttunen and I Vattulainen How would you integrate the equa tions of motion in dissipative particle dynamics simulations Com Phys Comm 153 407 2003 Igor Pasichnyk and Burkhard Diinweg Coulomb interactions via local dynamics A molecular
14. This cuboid in divided into bins_x bins in the x direction bins_y bins in the y direction and bins_z bins in the z direction such that the total number of bins is bins_x bins_y bins_z For each of these bins a stress tensor is calculated using the Irving Kirkwood method That is a given interaction contributes towards the stress tensor in a bin proportional to the fraction of the line connecting the two particles that is within the bin If the P3M and MMM1D electrostatic methods are used these interactions are not included in the local stress tensor The DH and RF methods in contrast are included Concerning bonded interactions only two body interactions FENE Harmonic are in cluded angular and dihedral are not For all electrostatic interactions only the real space part is included Care should be taken when using constraints of any kind since these are not accounted for in the local stress tensor calculations The command is implemented in parallel Output format variant 1 LocalStressTensor xbin y_bin z_bin y pressure_tensor specifying the local pressure tensor in each bin 8 2 Analyzing groups of particles molecules The following set of functions is designed to facilitate analysis of molecules Molecules are expected to be a group of particles comprising a contiguous range of particle IDs Each molecule is a set of consecutively numbered particles and all molecules are supposed to consist of the same
15. each particle has a likelihood of reacting in the vicinity of the surface distance is less than r as specified by the rate constant i e not according to Pey 1 exp nkAt with n the number of local catalysts To accomplish this each reactant is considered only once each time step by using the option react_once on The reaction command is set up such that the different properties may be influenced individually Only one type of particle can be assigned to each of these three reaction species and no particle type may be assigned to multiple species That is ESPResSo currently does not support particles of type 1 and 2 both to be reactants nor can particles of type 1 be a reactant as well as a catalyst Moreover only one of these reactions can be implemented in a single Tcl script If for instance there is a reaction involving particle types 1 2 and 4 there cannot be a second reaction involving particles of type 5 6 and 8 It is however possible to modify the reaction properties for a given set of types during the simulation Currently only strictly positive values of the catalytic conversion rate constant are allowed Setting the value to zero is equivalent to reaction off 86 e Variant 1 disables the reaction Note that at the moment there can only be one reaction in the simulation e Variant 2 returns the current reaction parameters In future versions of ESPResSo the capabilities of the CATALYTIC_REACTIONS
16. initarea 1 29 An initial guess for the area per lipid This guess is used to compute initial sphere dimensions based on the number of lipids This initial guess is then iteratively refined until all lipids can be fit uniformly on the sphere shuffle shuffle the topology prior to placing the lipids This is required for a random lipid distribution because otherwise the lipids will be placed on the sphere in the order they appear in the topology Creates a spherical vesicle by placing molecules in an ordered manner at uniform density on the surface of the sphere Molecules are assumed to have a uniform cross sectional area and closely matched though not identical lengths The radius of the vesicle will depend on the number of lipids and the area per lipid sphere_cap r arg half c arg initarea arg bondl arg shuffle r 10 0 The radius of the whole sphere where the cap is shaped half Create a half of sphere with the amount of molecules available c 0 0 0 0 0 0 The location of the center of the sphere relative to the center of the box initarea 1 29 An initial guess for the area per lipid This guess is used to compute initial sphere dimensions based on the number of lipids This initial guess is then iteratively refined until all lipids can be fit uniformly on the sphere shuffle shuffle the topology prior to placing the lipids This is required for a random lipid distribution because otherwise the lipids will be placed on the sph
17. particle distribution chains particles in the neighbourhood change_volume Tcl command pearl necklace structures configuration header pressure configure principal axis of the moment of in configure options 21 ertia 101 constraint Tcl commana 36 radial distribution function 110 copy_particles Tcl command radial distribution function g r 100 correlation Tcl command radius of gyration of a chain 108 Correlations 120 stress tensor 105 Coulomb interactions structure factor S q counterions Tcl command topologies 107 crosslink e van Hove autocorrelation function G r t 101 DAWAANR method Analysis in the Core Debye Hiickel potential analyze Tcl commana 97 diamond Tcl command Anisotropic interactions dielectric Tcl command 223 Dielectric interfaces dihedral interactions Dipolar interactions Directional Lennard Jones interaction 49 DLC method domain decomposition DPD DPD interaction dpd_gamma global variable dpd_r_cut global variable ELC method Electrostatic interactions end to end distance of a chain energies energy unit 9 features ADDITIONAL_CHECKS ADRESS ASYNC_ BARRIER BMHTF NACL BOND_ANGLE BOND_ANGLEDIST BOND_CONSTRAINT 198 BOND_ENDANGLEDIST BOND VIRTUAL BUCKINGHAM CATALYTIC REACTIONS CELL_DEBUG COLLISION DETECTION COMFIXED 198 COMFORCE COMM DEBUG CONSTRAINTS DIPOLES DPD 199 DPD_
18. pi the fluid velocity six floats Hse Uy yy Izz Ilyz II pi_neq the nonequilbrium part of the pressure ten sor components as above pop the 19 population check the order from the source code please boundary the flag indicating whether the node is a fluid node boundary 0 or a boundary node boundary 0 Does not support set Refer to the lbboundary command for this functionality Example The line puts lbnode 0 0 O print u prints the fluid velocity in node 0 0 0 to the screen The command set allows to change the density or fluid velocity in a single node Setting the other quantities can easily be implemented Example puts lbnode 0 0 O set u 0 01 0 0 12 4 Visualization Syntax lbfluid print vtk property filename Description The print parameter of the 1bfluid command is a feature to simplify visualization It allows for the export of the whole fluid field data into a file with name filename at once Currently supported values for the parameter property are boundary and velocity The additional option vtk enables export in the vtk format which is readable by visualization software such as paraview or mayaviZ Otherwise gnuplot readable data will be ex ported If you plan to use paraview for visualization note that also the particle positions can be exportet in the VTK format http www paraview org http code enthought com projects mayavi 151 12 5 Setting up boundary
19. script written using mbtools It is designed to run all of the examples provided but no more No doubt you will need to extend it for your own purposes 14 3 1 Variables used by main tcl main tcl expects the user to set various parameters in a parameters tcl file e g simplebilayer tcl Some of these parameters have defaults and generally don t need to be worried about except for specific cases In the following list variables that have no default and therefore must be set in the parameter file are noted with an asterisk e thermo Langevin The type of thermostat to be used Set to DPD for a dpd thermostat Any other value gives a langevin e dpd_gamma Required if you set the thermo to DPD e dpd_r_cut Required if you set the thermo to DPD e warmup_temp systemtemp The temperature at which the warmup is performed The default behaviour is to use the system temperature e warmsteps 100 Number of integrate steps per warmup cycle e warmtimes 20 Number of calls to integrate over which the warmup occurs e free_ warmsteps 0 Warmup steps to be used for the warmup that occurs after particles are freed of any temporary constraints e free_warmtimes 0 Warmup times to be used for the warmup that occurs after particles are freed of any temporary constraints 160 npt off Whether to use the constant pressure barostat p_ext The pressure you want to simulate at Required if npt is set to on piston_mass box mass Required if npt
20. 104 range_start_y range_start_z range_x range_y range_z bins_x bins_y bins _z analyze set chains chain_start n_chains chain_length 105 analyze set topo_part_sync analyze set analyze re lt re gt chain_start n_chains chain_length 106 analyze rg lt rg gt chain_start n_chains chain_length 107 analyze rh lt rh gt chain_start n_chains chain_length 108 analyze internal_dist lt internal_dist gt chain_start n_chains chain_1d108 analyze bond_dist lt bond_dist gt index index 108 chain_start n_chains chain_length analyze bond_1 lt bond_1 gt chain_start n_chains chain_length 109 analyze formfactor lt formfactor gt qmin qmaz qbins 109 chain_start n_chains chain length analyze rdfchain rmin rmaz rbins chain_start n_chains chain_length 110 analyze lt gi gt lt g2 gt lt g3 gt chain_start n_chains chain_length 110 analyze g123 init chain_start n_chains chain_length 19 ww analyze append analyze remove index analyze replace index analyze push size analyze configs config analyze configs analyze stored uwerr data nrep col s_tau plot uwerr data nrep f s_tau f_args plot observable new name parameters observable id print formatted observable id delete observable new needs_profile_specs other_parameters minx minz maxx mazz miny miny maxy mary minz minz maxz mazz xbins zbins ybins ybins zbins zbins
21. 2 LB as a thermostat Syntax thermostat 1b T Required features LB LB_GPU Description The LBM implementation in ESPResSo uses Duenweg s point coupling method to couple MD particles the LB fluid This coupling consists in a frictional force and a random force gt F y T Er The frictional force tends to decrease the relative velocity between the fluid and the particle whereas the random forces are chosen so large that the average kinetic energy per particle corresponds to the given temperature according to a fluctuation dissipation theorem No other thermostatting mechanism is necessary then Please any of these off before starting the LB thermostatting mechanism The LBM implementation provides a fully thermalized LB fluid i e all nonconserved modes including the pressure tensor fluctuate correctly according to the given temper ature and the relaxation parameters All fluctuations can be switched off by setting the temperature to 0 12 3 Reading and setting single lattice nodes Syntax lbnode z y z print set args Required features LB Description The 1bnode command allows to inspect print and modify set single LB nodes Note that the indexing in every direction starts with 0 For both commands you have to specify what quantity should be printed or modified Print allows the following arguments 150 rho the density scalar u the fluid velocity three floats uz Uy uz
22. 2 imd positions unfolded fold_chains 3 imd listen seconds 4 imd disconnect Description In your simulation the IMD connection is setup up using variant 1 where port is an arbitrary port number which usually has to be between 1024 and 65000 By default ESPResSo will try to open port 10000 but the port may be in use already by another ESPResSo simulation In that case it is a good idea to just try another port While the simulation is running variant 2 can be used to transfer the current co ordinates to VMD if it is connected If not nothing happens and the command just consumes a small amount of CPU time Note that before you can transfer coordinates to VMD VMD needs to be aware of the structure of the system For that you first need to load a corresponding structure file PSF or VSF into VMD Also note that the command prepare_vmd_connection see section 10 6 3 on the next page can be used to automatically set up the VMD connection and transfer the structure file By specifying unfolded the unfolded coordinates of the particles will transferred while fold_chains will fold chains according to their centers of mass and retains bond ing connectivity Note that this requires the chain structure to be specified first using the analyze command 136 Variant 3 can be used to let the simulation wait for seconds seconds or until IMD has connected before the script is continued This is normally only useful in
23. 3 checks for a block with tag block and then again executes the corresponding blockfile_read_auto_tag if it exists In the contrary that means that for a new blocktype you will normally implement two procedures blockfile_write_tag channel write tag arg which writes the block including the header and enclosing braces and blockfile_read_auto_tag channel read auto which reads the block data and the closing brace The parameters write read tag and auto are regular parameters which will always have the specified value They occur just for technical reasons In a nutshell The blockfile command is provided for saving and restoring the current state of ESPResSo e g for creating and using checkpoints Hence you can transfer all accessible information from files to ESPResSo and vice versa set out open lgzip c gt checkpoint block gz w blockfile out write variable all blockfile out write interactions blockfile out write random blockfile out write bitrandom blockfile out write particles id pos type q v f all blockfile out write bonds all blockfile out write configs close out This example writes all global variables all interactions the full current state of the random number generator all information i e id position type number charge velocity forces bonds of all particles and all stored particle configurations to the file checkpoint block gz which is compressed on the fly If you want to be able
24. 8 tetra functional nodes e monomers_per_chain Sets the number of chain monomers between the functional nodes e counterions Nc Adds Ncr counterions to the system 33 e charges Valnode V lnonomer Valci Sets the charge of the nodes to Valnode the charge of the connecting monomers to valmonomer and the charge of the counte rions to valci e distance dehargea Specifies the distance between charged monomers along the interconnecting chains If dehargea gt 1 the remaining chain monomers are un charged e nonet Do not create bonds between the chains 4 2 5 icosaeder Setting up an icosaeder Syntax icosaeder a monomers_per_chain counterions Nc 1 1 charges Valmonomers Valci distance decharged Required features ELECTROSTATICS Description Creates a modified icosaeder to model a fullerene or soccer ball The edges are modeled by polymer chains connected at the corners of the icosaeder For inter particle bonds interaction O is taken which must be a two particle bond Two particle types are used for the pentagons and the interconnecting links For an example see figure 4 2 Figure 4 2 Icosaeder with monomers_per_chain 15 Arguments ea Length of the links Defines the size of the icosaeder e monomers_per_chain Specifies the number of chain monomers along one edge e counterions Nc Specifies the number of counterions to be placed into the sys tem 34 e charges Valmonomers Valci
25. DPD either via a global thermo stat or via a thermostat and a special DPD interaction between particle types The latter allows the user to specify friction coefficients on a per interaction basis Thermostat DPD Syntax thermostat dpd temperature gamma r_cut WF wf tgamma tr_cut TWF twf Required features DPD or TRANS_DPD Description ESPResSo s standard DPD thermostat implements the thermostat exactly as described in 50 We use the standard Velocity Verlet integration scheme e g DPD only influences the calculation of the forces No special measures have been taken to self consistently determine the velocities and the dissipative forces as it is for example described in 87 DPD adds a velocity dependent dissipative force and a random force to the usual con servative pair forces e g Lennard Jones 77 The dissipative force is calculated by D D oo Fig Cw rig Pig Vig Pig The random force by AR R a Fij FU rij OijPij where O 0 5 0 5 is a uniformly distributed random number The connection of o and is given by the dissipation fluctuation theorem ow ri Cw rij kpT The parameters gamma r_cut define the strength of the friction and the cutoff radius According to the optional parameter WF can be set to 0 or 1 default is 0 of the thermostat command the functions w and w are chosen in the following way Tus r_cut EE ee r WED Pr wry f 1 _WF 1 For rij gt
26. For a typical queueing system a starting routine could look like this master for h in HOSTS do if master then ssh h cd run pt_test tcl master h else ssh h cd run pt_test tcl connect host fi done where pt_test tcl passes the connect option on to parallel_tempering main Sharing data Syntax parallel_tempering set_shareddata data Description can be used at any time by the master process to specify additional data that is available on all processes of the parallel_tempering simulation The data is accessible from all processes as parallel_tempering shareddata 94 7 7 Metadynamics Syntax 1 metadynamics 2 metadynamics set off 3 metadynamics set distance pid pido dmin dmax Oneight Dwidth fbound dbins 4 metadynamics set relative_z pid pidg 2min 2max Dheight width found 2bins 5 metadynamics print_stat current_coord 6 metadynamics print_stat coord_values 7 metadynamics print_stat profile 8 metadynamics print_stat force 9 metadynamics load_stat profile_list force_list Required features METADYNAMICS Description Performs metadynamics sampling Metadynamics is an efficient scheme to calculate the potential of mean force of a system as a function of a given reaction coordinate from a canonical simulation The user first chooses a reaction coordinate e g distance between two particles pid and pida As the system samples values along this reaction coordina
27. For real simulations the simulation scripts can extend over thousands of lines of code and contain automated adaption of parameters or online analysis up to automatic generation of data plots Parameters can be changed arbitrarily during the simulation process as needed for e g 18 simulated annealing The possibility to perform non standard simulations without the need of modifications to the simulation core was one of the main reasons why we decided to use a script language for controlling the simulation core 2 4 tutorial tcl In the directory samples of the es sources you will find a well documented simulation script tutorial tcl which takes you step by step through a slightly more complicated simulation of a polyelectrolyte system The basic structure of the script is however the same as in the previous example and probably the same as the structure of most ESPResSo simulation scripts Initially some parameters and global variables are set the interactions are initialized and particles are added For this the script makes use of the polymer command which provides a faster way to set up chain molecules The actual simulation falls apart again into two loops the warmup loop with increasing force capping and the final simulation loop Note that the electrostatic interaction is only activated after equilibrating the excluded volume interactions which speeds up the warmup phase However depending on the problem this splitted
28. How often to print out configurations e umdcommands a list of additional lines of commands to be written to the vmd_ animation script file 14 4 Analysis The analysis package is designed to help organise the many possible analysis routines that can be performed during a simulation This documentation describes the basic user interface commands and then all of the possible analysis routines Instructions on how to add a new analysis routine are given at the end of this section 14 4 1 Basic commands The following commands comprise the user interface to the analysis package At the start of a simulation all of the analysis that is to be performed is specified using the setup_analysis command Each time you want the analysis performed a call to do_analysis should be made One can then call print_averages to write results to file mbtools analysis setup_analysis outputdir arg suffix arg iotype arg g arg str arg e commands A tcl list where each element of the list is a string specifying the name and complete argument list for a particular analysis to be carried out e outputdir The directory where analysis output files will be created e suffix tmp Suffix that will be appended to standard file names for analysis output e iotype a The iotype that will be used when opening files for analysis For an explanation of the different iotypes see the documentation for the standard tcl command open e y 8 Number of grid
29. In ESPResSo the LB scheme and the MD scheme are not synchronized In one LB time step typically several MD steps are performed This allows to speed up the simulations and is adjusted with the parameter tau the LB timestep The parameters dens and visc set up the density and kinematic viscosity of the LB fluid in usual MD units 148 Internally the LB implementation works with a different set of units all lengths are expressed in agrid all times in tau and so on Therefore changing agrid and tau might change the behaviour of the LB fluid e g at boundaries due to characteristics of the LBM itself It should also be noted that the LB nodes are located at 0 5 1 5 2 5 etc in terms of agrid This has important implications for the location of hydrodynamic boundaries which are generally considered to be halfway between two nodes to first order Currently it is not possible to precisely give a parameter set where reliable results are expected but we are currently performing a study on that Therefore the LBM should not be used as a black box but only after a careful check of all parameters that were applied The parameter ext_force allows to apply an external body force density that is homogeneous over the fluid It is again to be given in MD units The parameter bulk_viscosity allows to tune the bulk viscosity of the fluid and is given in MD units In the limit of low Mach often also low Reynolds number the results should be inde pen
30. Set the charges of the monomers to valmonomers and the charges of the counterions to valci e distance dehargea Specifies the distance between two charged monomer along the edge If dchargea gt 1 the remaining monomers are uncharged 4 2 6 crosslink Cross linking polymers Syntax crosslink num_polymer monomers_per_chain start pid catch Teatch distLink link_dist distChain chain dist FENE bondid trials trymax Description Attempts to end crosslink the current configuration of num_polymer equally long poly mers with monomers_per_chain monomers each returning how many ends are success fully connected Arguments e start pid pid specifies the first monomer of the chains to be linked It has to be specified if the polymers do not start at id 0 e catch ratch Set the radius around each monomer which is searched for possible new monomers to connect to Teatch defaults to 1 9 e distLink link_dist The minimal distance of two interconnecting links It de faults to 2 e distChain chain_dist The minimal distance for an interconnection along the same chain It defaults to 0 If set to monomers_per_chain no interchain con nections are created e FENE bondid Sets the bond type for the connections to bondid e trials trymax If not specified trymax defaults to 30000 4 2 7 copy_particles copying a set of particles Syntax copy_particles set dl id2 range from to shift sx sy s z Descript
31. Xx a a The parameters x k 1 Kk 1 and x ag z 1 o 1 are responsi ble for the degree of anisotropy of the molecular properties ky is the molecular elonga tion and kz is the ratio of the potential well depths for the side by side and end to end 51 configurations The exponents u and y are adjustable parameters of the potential Sev eral Gay Berne parametrizations exist the original one being ky 3 kg 5 u 2 and v 1 5 3 Bonded interactions Syntax inter bondid interaction parameters Description Bonded interactions are identified by their bonded interaction type identificator bondid which is a non negative integer The inter bondid command is used to specify the type and parameters of a bonded interaction which applies to all par ticles connected explicitely by this bond using the part command see section 4 1 on page 26 Therefore defining a bond between two particles always involves two steps defining the interaction and applying it Assuming that two particles with ids 42 and 43 already exist one can create e g a FENE bond between them using inter 1 fene 10 0 2 0 part 42 bond 1 43 If a FENE bond with the same interaction parameters is required between several particles e g in a simple chain molecule one can use the sampe type id inter 1 fene 10 0 2 0 part 42 bond 1 43 part 43 bond 1 44 Bonds can have more than just two bond partners For the inter command tha
32. a product de composition or similar tricks MMM1D has to be implemented as a simple NxN loop However the formulas can be evaluated efficiently so that MMMID can still be used reasonably for up to 400 particles on a single processor 3 E 4 ELC The ELC method differs from the other MMM algorithms in that it is not an algorithm for the calculation of the electrostatic interaction but rather represents a correction term which allows to use any method for threedimensionally periodic systems with spherical summation order for twodimensional periodicity The basic idea is to expand the two dimensional slab system of height h in the non periodic z coordinate to a system with periodicity in all three dimensions with a period of A gt h which leaves an empty gap of height 6 A h above the particles in the simulation box Since the electrostatic potential is only finite if the total system is charge neutral the additional image layers those layers above or below the original slab system are charge neutral too Now let us consider the n th image layer which has an offset of nA to the original layer If nA is large enough each particle of charge q_j at position 4 yj Zj n z and its replicas in the xy plane can be viewed as constituting a homogeneous charged sheet of charge density oj Oe The potential of such a charged sheet at distance z is 2705 2 Now we consider the contribution from a pair of image layers located at nA n
33. and the ratio between radii cylinder c arg initarea arg bondl arg shuffle e c 0 0 0 0 0 0 e initarea 1 29 e shuffle shuffle the topology prior to placing the lipids Creates a cylinder which spans the box along one dimension by placing molecules uni formly on its surface Works in a similar way to the sphere routine random exclude arg inside arg shuffle bondl arg e exclude arg an exclusion zone definition suitable for passing to mbtools utils isoutside e inside arg an inclusion zone definition suitable for passing to mbtools utils isoutside e shuffle shuffle the topology prior to placing the lipids 170 Places molecules randomly in space with a sortof random orientation vector If an exclusion zone is defined then no molecules will be placed such that their positions are within the zone If an inclusion zone if defined then no molecules will be place outside this zone For instance set geometry geometry random exclude sphere 0 0 0 0 0 0 4 0 inside cuboid 0 0 0 0 0 0 15 0 15 0 15 0 will randomly place molecules within the volume between a sphere with a radius of 4 0 and a cuboid with dimension 15 0 x 15 0 x 15 0 at the origin readfile ignore arg f arg t arg e ignore arg particle properties to be ignored during the file read f arg The file containing the configuration to be used for setup Must be an espresso blockfile with box length part
34. artificially cap the forces which often allows to equilibrate the system much faster See the subsection for details 43 5 1 1 Tabulated interaction Syntax inter typel type2 tabulated filename Required features TABULATED Description This defines an interaction between particles of the types type1 and type2 according to an arbitrary tabulated pair potential filename specifies a file which contains the tabulated forces and energies as a function of the separation distance The tabulated potential allows capping the force using inter forcecap see section 5 9 5 At present the required file format is simply an ordered list separated by whitespace The data reader first looks for a character and begins reading from that point in the file Anything before the will be ignored The first three parameters after the specify the number of data points Npoints and the minimal and maximal tabulated separation distances Tmin and Tmax The number of data points obviously should be an integer the two other can be arbitrary positive doubles Take care when choosing the number of points since a copy of each lookup table is kept on each node and must be referenced very frequently The maximal tabulated separation distance also acts as the effective cutoff value for the potential The remaining data in the file should consist of n data triples r F r and V r r gives the particle separation V r specifies the interaction potential and F r
35. averaged over all stored configurations see section 8 3 on page 112 This function assumes linear chain topology and does not check if the bonds really exist Output format mean stddev max min Form factor Syntax analyze formfactor lt formfactor gt gmin qmax qbins chain_start n_chains chain length Description Computes the spherically averaged form factor of a single chain which is defined by chain_length Y sin qri 8 6 pga Ww 1 Pg chain_length of a single chain averaged over all chains for gbin 1 logarithmically spaced q vectors qmin qmaz where qmin gt 0 and qmaz gt qmin If lt formfactor gt is used the form factor will be averaged over all stored configurations see section 8 3 on page 112 Output format 4 Sla with q qmin qmaz Chain radial distribution function Syntax analyze rdfchain rmin rmazx rbins chain_ start n_chains chain_length 110 Description Returns three radial distribution functions rdf for the chains The first rdf is calculated for monomers belonging to different chains the second rdf is for the centers of mass of the chains and the third one is the distribution of the closest distances between the chains i e the shortest monomer monomer distances The distance range is given by rmin and rmaz and it is divided into rbins equidistant bins Output format 1 r rdf1 r rdf2 r rdf3 r y Mean square displacement o
36. blockfile_write_tag is called with all of the commands arguments This function should then write the data Example set file open data dat w blockfile file write start mydata puts file This is my data blockfile file write end will write mydata This is my data to the file data dat 10 2 7 Reading blocks Syntax 1 blockfile channel read start 2 blockfile channel read toend 3 blockfile channel read particles interactions bonds variable seed random bitrandom configs 4 blockfile channel read auto Description Variants 1 and 2 are the low level block reading commands Variant 1 reads the start part of a block and returns the block title while variant 2 reads the block data and returns it Variants 3 and 4 read whole blocks Variant 3 reads the beginning of one block checks wether it contains data of the given type and reads it Variant 4 reads in one block and does the following 1 if a procedure blockfile_read_auto_tag exists this procedure takes over tag is the first expression in the block For most block types at least all mentioned above i e particles interactions bonds seed random bitrandom configs and variable the corresponding procedure will overwrite the current information with the information from the block 130 2 if the procedure does not exist it returns usertag rest_of_block 3 if the file is at the end it returns eof Variant
37. bond without associated potential or force It can used to specify topologies and for some analysis that rely on bonds or e g for bonds that should be displayed in VMD 54 5 4 Object in fluid interactions Please cite BIBTEX key cimrak in file doc ug citations bib when using the interactions in this section in order to simulate extended objects embedded in a LB fluid The following interactions are implemented in order to mimic the mechanics of elastic or rigid objects immersed in the LB fluid flow Their mathematical formulations have been taken from 18 Details on how the bonds with fluid structure interactions can be used and automated are described in section 5 4 1 Stretching force Syntax inter bondid stretching_force L p ks Description This type of interaction is available for closed 3D immersed objects as well as for 2D sheet flowing in the 3D flow For each edge of the mesh Lap is the current distance between point A and point B L9 p is the distance between these points in the relaxed state that is if the current edge has the length exactly En pg then no forces are added ALag is the deviation from the relaxed state that is ALAB Lap En pg The stretching force between A and B is computed using ALAB 0 Lab F A B kkl Aan HAB 5 21 Here nag is the unit vector pointing from A to B ks is the stretching constant AAB LaB L g and is a nonlinear function that resembles n
38. by the user You can find some examples in the samples folder of the source code directory If you want to run in parallel you should have compiled ESPResSo to use MPI and need to tell MPI to run ESPResSo in parallel The actual invocation is implementation dependent but in many cases such as OpenMPI you can use mpirun n n_nodes Espresso script where n_nodes is the number of prcessors to be used 2 2 Running ESPResSo ESPResSo is implemented as an extension to the Tcl scripting language This means that you need to write a script for any task you want to perform with ESPResSo To learn about the Tcl script language and especially the ESPResSo extensions this chapter offers two tutorial scripts The first will guide you step by step through creating your first simulation script while the second script is a well documented example simulation script Since the latter is slightly more complex and uses more advanced features of ESPResSo we recommend to work through both scripts in the presented order If you want to learn about the Tcl language in greater detail there is an excellent tutorial 2 3 Creating the first simulation script This section introduces some of the features of ESPResSo by constructing step by step a simulation script for a simple salt crystal We cannot give a full Tcl tutorial here however most of the constructs should be self explanatory We also assume that the reader is familiar with the basic concepts of a M
39. by typing make dg in the build directory Afterwards it can be found in the subdirectory of the build directory doc dg html index html A recent version of this guide can also be found on the ESPResSo homepage espressomd org 16 4 User s guide If while reading this user guide you notice any mistakes or badly if at all described features or commands you are very welcome to contribute to the guide and have others benefit from your knowledge For this you should also checkout the development version as described on the home page As the user guide like all ESPResSo code is always in flow and changes are made regularly there are already many paragraphs marked with a todo box To turn on these boxes edit the main file doc ug ug tex and adapt the inclusion of the ATEX package todonotes You can then build the user guide by typing make ug 185 A ESPResSo quick reference part pid pos x y z type typeid v vz vy vz fe fy fz 26 bond bondid pid q charge quat q1 q2 q3 q4 omega x y z torque z y z rinertia x y z un fix x y z ext_force z y 2 ext_torque z y z exclude pid2 exclude delete pid2 mass mass dipm moment dip dx dy dz virtual v vs_relative pid distance vs_auto_relate_to pid temp T gamma g rotation rot Required features ELECTROSTATICS ROTATION EXTERNAL_FORCES EXCLUSION 9mass DIPOLES VIRTUAL_SITES_COM VIRTUAL_SITES_RELATIVE
40. chain ebond_length Sets the initial distance between two adjacent monomers The dis tance during the course of the simulation depends on the applied potentials For fixed bond length please refer to the Rattle Shake algorithm 2 The algorithm is based on Verlet algorithm and satisfy internal constraints for molecular models with internal constrains using Lagrange multipliers start pid Sets the particle number of the start monomer to be used with the part command This defaults to 0 e pos x y z Sets the position of the first monomer in the chain to zx y z defaults to a randomly chosen value mode RW PSAW SAW shield trymax Selects the setup mode 30 RW Random walk The monomers are randomly placed by a random walk with a steps size of bond ength PSAW Pruned self avoiding walk The position of a monomer is randomly cho sen in a distance of bond _length to the previous monomer If the position is closer to another particle than shield the attempt is repeated up to trymax times Note that this is not a real self avoiding random walk as the particle distribution is not the same If you want a real self avoiding walk use the SAW mode However PSAW is several orders of magnitude faster than SAW especially for long chains SAW Self avoiding random walk The positions of the monomers are chosen as in the plain random walk However if this results in a chain that has a monomer that is closer to another pa
41. corresponding particle The membrane is thus discretized into mesh_nnode particles with IDs starting from 0 to mesh_nnode 1 The IDs are assigned in the same order as in the mesh nodes dat file The mesh triangles dat contains mesh_ntriangle lines with three nonnegative in tegers separated by blank space Each line represents one triangle in the triangulation For algorithmic purposes it is crucial to have defined a correct orientation of the trian gle The orientation is defined using the normal vector associated with the triangle The important rule is that the normal vector of the triangle must point inside the immersed object As an example let us have one line in the file mesh triangles dat with numbers 4 0 and 7 This means that particles with IDs 4 0 and 7 form one triangular face of the triangulation The orientation is defined as follows create two vectors v and v2 such that v is pointing from particle 4 to particle 0 and va is pointing from particle 4 to particle 7 Be carefull the order of vectors and particles matters The normal vector n is computed as a vector product v x v2 The direction of n can be determined by the rule of right hand the thumb points in the v direction the index finger in the vz direction and the middle finger in the n direction Following this principle all the lines in the mesh triangles dat files must be such that the normal vectors of the corresponding triangles must point inside the immersed ob
42. directory itself but to start configure from another directory see section 3 1 1 on page 21 Furthermore many features of ESPResSo can be selectively turned on or off in the local configuration header see section 3 4 on page 24 before starting the compilation with make The shell script configure prepares the source code for compilation It will determine how to use and where to find the different libraries and tools required by the compilation process and it will test what compiler flags are to be used The script will find out most of these things automatically If something is missing it will complain and give hints on how to solve the problem The configuration process can be controlled with the help of a number of options that are explained in section 3 1 on page 20 The command make will compile the source code Depending on the options passed to the program make can also be used for a number of other things e It can install and uninstall the program to some other directories However nor mally it is not necessary to actually install ESPResSo to run it e It can test ESPResSo for correctness e It can build the documentation The details of the usage of make are described in section 3 2 on page 22 When these steps have successfully completed ESPResSo can be started with the command see section 3 3 on page 23 Espresso script 12 where script is a Tcl script that tells ESPResSo what to do and has to be written
43. e g be seen when using the sample script poisseuille tcl with a high viscosity The bounce back boundary conditions allow to set velocity at a boundary to a nonzero value This allows to create shear flow and boundaries moving relative to each other This could be a fixed sphere in a channel moving at a finite speed corresponding to the galilei transform of a moving sphere in a fixed channel The velocity boundary conditions are implemented according to eq 12 58 Using this implementation as a blueprint for the boundary treatment an implementation of the Ladd Coupling should be relatively straightforward Syntax lbboundary force Nboundary Required features LB_BOUNDARIES 152 Description This variant prints out the force on boundary number n_boundary 12 6 Choosing between the GPU and CPU implementations Syntax 1 lbfluid cpu 2 lbfluid gpu Required features lig 2 LB_GPU Description A very recent development is an implementation of the LBM for NVIDIA GPUs using the CUDA framework On CUDA supporting machines this can be activated by config uring with configure with cuda path to cuda and activating the feature LB_GPU Within the ESPResSo Tcl script the 1bfluid command can be used to choose between the CPU and GPU implementations of the Lattice Boltzmann algorithm for further information on CUDA support see section 6 6 Variant 1 is the default and turns on the standard CPU implementation of the Lattice Bol
44. feature may be generalized to handle multiple reactant catalyzer and product types as well as more general reaction schemes Other changes may involve merging the current implementation with the COLLISION_DETECTION feature 6 9 Galilei Transform and Particle Velocity Manipulation The following commands may be useful in effecting the velocity of the system 6 9 1 Particle motion and rotation Syntax kill_particle_motion rotation Required features ROTATION Description This command halts all particles in the current simulation setting their velocities to zero as well as their angular momentum if the option rotation is specified and the feature ROTATION has been compiled in 6 9 2 Forces and torques acting on the particles Syntax kill_particle_forces torques Required features ROTATION Description This command sets all forces on the particles to zero as well as all torques if the option torque is specified and the feature ROTATION has been compiled in 6 9 3 The centre of mass of the system Syntax system_CMS Description Returns the center of mass of the whole system It currently does not factor in the density fluctuations of the Lattice Boltzman fluid 87 6 9 4 The centre of mass velocity Syntax system_CMS_velocity Description Returns the velocity of the center of mass of the whole system 6 9 5 The Galilei transform Syntax galilei_transform Description Substracts the ve
45. following meaning epsilon epsilon The dielectric constant of the surrounding medium metallic i e in finity or some finite positive number Defaults to metallic n_interpol n nterpol Number of interpolation points for the charge assignment func tion When this is set to 0 interpolation is turned off and the function is computed directly Defaults to 32768 mesh_off mesh ff Offset of the first mesh point from the lower left corner of the simu lation box in units of the mesh constant Defaults to 0 5 0 5 0 5 5 7 2 Debye Htckel potential Syntax inter coulomb lg dh K Teut Required features ELECTROSTATICS Description Defines the electrostatic potential by UC PH lgkgT for r lt Teut 5 34 q uqerp r D The Debye H ckel potential is an approximate method for calculating electrostatic interactions but technically it is treated as other short ranged non bonding potentials For r gt freut it is set to zero which introduces a step in energy Therefore it introduces fluctuations in energy For 0 this corresponds to the plain coulomb potential 63 5 7 3 MMM2D Please cite 5 BIBTEX key mmm2d in file doc ug citations bib when using MMM2D and 55 BIBTEX key icmmm2d in file doc ug citations bib when using dielectric interfaces Syntax inter coulomb lg mmm2d mazimal_pairwise_error fixed_far_cutoff dielectric m p dielectric contrasts A Ap Required features ELEC
46. for non bonded tabulated potentials see section 5 1 1 The potential is calculated as follows e Variant 1 is a two body interaction depending on the distance of two particles The force acts in the direction of the connecting vector between the particles The bond breaks above the tabulated range but for distances smaller than the tabulated range a linear extrapolation based on the first two tabulated force values is used e Variant 2 is a three body angle interaction similar to the angle potential see section B 5 It is assumed that the potential is tabulated for all angles between 0 and 7 where 0 corresponds to a stretched polymer and just as for the tabu lated pair potential the forces are scaled with the inverse length of the connecting vectors The force on particles p and p3 in the notation of section 5 5 acts per pendicular to the connecting vector between the particle and the center particle p2 in the plane defined by the three particles The force on the center particle pa balances the other two forces e Variant 3 tabulates a torsional dihedral angle potential see section 5 6 It is assumed that the potential is tabulated for all angles between 0 and 27 This potential is not tested yet Use on own risk and please report your findings and eventually necessary fixes 5 3 6 Virtual bonds Syntax inter bondid virtual_bond Description This creates a virtual bond type with identificator bondid i e a pair
47. from the origin in units of the normal vector so that the product of d and n is a point on the surface Therefore negative distances are quite common The resulting surface in variant 2 is a sphere with center cy c and radius rad The direction determines the force direction 1 or inside for inward and 1 or outside for outward The resulting surface in variant 3 is a cylinder with center cz cy cz and radius rad The length parameter is half of the cylinder length The axis is a vector along the cylinder axis which is normalized in the program The direction is defined the same way as for the spherical constraint The resulting surface in variant 4 is a rhomboid defined by one corner located at Px Py pz and three adjacent edges defined by the three vectors connecting the corner p with it s three neighboring corners a az ay az b bz by bz and c cz Cy Cz The resulting surface in variant 5 is n spheres of radius r along each dimension connected by cylinders of radius re The spheres have simple cubic symmetry The spheres are distributed evenly by dividing the box by n Dimension of the maze can be controlled by d 0 for one dimensional 1 for two dimensional and 2 for three dimensional maze Variant 6 sets up a cylindrical pore similar to variant 3 with a center cz Cy c and radius rad The length parameter is half of the cylinder length The azis is a vector 37 along the cylinder axis which is normalized
48. in the program The argument radius rad can be replaced by the argument radii radi rad2 to obtain a pore with a conical shape and corresponding opening radii The first radius is in the direction opposite to the axis vector Variant 7 specifies an electrostatic interaction between the charged particles in the system to an infinitely long rod with a line charge of lambda which is alinge along the z axis and centered at Cz and Cy Variant 8 specifies the electrostatic interactinos between the charged particles in the system and an inifinitely large plate in the x y plane at height h The plate carries a charge density of sigma Variant 9 specifies the dipolar coupling of particles with a dipolar moment to an external field fz fy fz Variant 10 creates an infinite plane at a fixed position For non initializing a direc tion of the constraint values of the positions have to be negative For the tunable slip boundary interactions you have to set two constraints Example To create an infinite plane in z direction at z 20 0 of type id 1 use constraint plane cell 10 10 20 type 1 4 3 1 Deleting a constraint Syntax constraint delete num Description This command will delete constraints If num is specified only this constraint will deleted otherwise all constraints will be removed from the system 4 3 2 Getting the force on a constraint Syntax constraint force n Description Returns the force acting on the nth c
49. inclusion of hydrodynamic interac tions into the simulation The implementation of boundary conditions for the LBM is a difficult task a lot of research is still being conducted on this topic The focus of the ESPResSo implementation of the LBM is of course the coupling to MD and therefore available geometries and boundary conditions are somewhat limited in comparison to pure codes Here we restrict the documentation to the interface For a more detailed description of the method please refer to the literature 12 1 Setting up a LB fluid Please cite 7 BIBTEX key espresso2 in file doc ug citations bib if you use the LB fluid and BIBTEX key 1bgpu in file doc ug citations bib if you use the GPU implementation Syntax l or 2 JE ets dens density bulk_visc bulk _viscosity l or 2 lbfluid gpu agrid agrid tau lb timestep ext_force f fy f friction gamma l or 2 visc viscosity l or gamma_odd gamma_odd l or 2 lip 218 GPU gamma_even gamma_even Required features Description The 1bfluid command initializes the fluid with a given set of parameters It is also possible to change parameters on the fly but this will only rarely be done in practice Before being able to use the LBM it is necessary to set up a box of a desired size The parameter agrid is used to set the lattice constant of the fluid so the size of the box in every direction must be a multiple of agrid
50. its top layer s to the next processor dealing with the layers above Simultaneously the top processor starts with the calculation of the alhs 0 0 and sends them down After the communicated has been completed every processor can use the E ae and the eae contributions for its particles In pseudo code the far formula algorithm looks like to calculate the force rsp energy 1 for each layer s 1 5 s c s c 0 ms a E b for each particle j in layer s E s c s c i calculate ES 1 A ue ace ere l s c s c s c s c 9 al el _ a bs es e 3 for each layer s 4 S 213 a a 0 8 c 5 c ie s s 2 all s c s c s e 8 c 4 S 2 S 5 for each layer s S 3 1 l s c s c l s c s c s c s c a Elis es e glls es o y sese 6 for each layer s 1 S a for each particle j in layer s o E aid See E i calculate particle interaction from For further details see Arnold and Holm 5 4 Arnold et al 6 E 2 1 Dielectric contrast A dielectric contrast at the lower and or upper simulation box boundary can be included comparatively easy by using image charges Apart from the images of the lowest and topmost layer the image charges are far enough to be treated by the far formula and can be included as starting points in the calculation of the terms The remaining particles from the lowest and topmost layer are tre
51. more details on the format refer to the homepage of the formad Creating files in these formats from within ESPResSo is supported by the commands writevsf and writevcf that write a structure respectively a coordinate block to the given Tcl channel To create a VTF file first use writevsf at the beginning of the simulation and then writevcf after each timestep to generate a trajectory of the whole simulation The structure definitions in the VTF VSF formats are incremental i e a user can easily add further structure lines to the VTF VSF file after a structure block has been written to specify further particle properties for visualization Note that the ids of the particles in ESPResSo and VMD may differ VMD requires the particle ids to be enumerated continuously without any holes while this is not required in ESPResSo When using writevsf and writevcf the ESPResSo particle ids are automatically translated into VMD particle ids The function vtfpid allows the user to get the VMD particle id for a given ESPResSo particle id Also note that these formats can not be used to write trajectories where the number of particles or their types varies between the timesteps This is a restriction of VMD itself not of the format 10 3 1 writevsf Writing the topology Syntax writevs channelld short verbose radius radii auto typedesc typedesc ignore_charges Description Writes a structure block describing the system s struc
52. not be discussed here but the penalty is at least 10ns resulting in an effective memory transfer rate of only 800MB s To remedy this modern processors have a small amount of low latency memory directly attached to the processor the cache The processor cache is organised in different levels The level 1 L1 cache is built directly into the processor core has no latency and delivers the data immediately on demand but has only a small size of around 128KB This is important since modern processors can issue several simple operations such as additions simultaneously The L2 cache is larger typically around 1MB but is located outside the processor core and delivers data at the processor clock rate or some fraction of it In a typical implementation of the link cell scheme the order of the particles is fairly random determined e g by the order in which the particles are set up or have been communicated across the processor boundaries The force loop therefore accesses the particle array in arbitrary order resulting in a lot of unfavourable page misses In the memory organisation of ESPResSo the particles are accessed in a virtually linear order Because the force calculation goes through the cells in a linear fashion all accesses to a single cell occur close in time for the force calculation of the cell itself as well as for its neighbours Using the domain decomposition cell scheme two cell layers have to be kept in the processor cache For 1
53. number of particles Some functions in this group require that the particles constituting a molecule are connected into linear chains particle n is connected to n 1 and so on while others are applicable to molecules of whatever topology The analyze set command defines the structure of the current system to be used with some of the analysis functions Syntax 1 analyze set chains chain_start n_chains chain length 2 analyze set topo_part_sync 3 analyze set Description Variant 1 defines a set of n_chains chains of equal length chain_length which start with the particle with particle number chain_start and are consecutively numbered i e the last particle in that topology has number chain_start n_chains x chain_length 1 Variant 2 synchronizes topology and particle data assigning mol_id values to parti cles 107 Variant 3 will return the chains currently stored 8 2 1 Chains All analysis functions in this section require the topology of the chains to be set correctly The topology can be provided upon calling This re sets the structure info permanently i e it is only required once End to end distance Syntax analyze re lt re gt chain_start n_chains chain length Description Returns the quadratic end to end distance and its root averaged over all chains If lt re gt is used the distance is averaged over all stored configurations see section age 17 Output format re error_of_
54. of mean force of the end to end distance of a poly mer Variant 4 sets a metadynamics scheme with the reaction coordinate relative_z relative height i e z coordinate of particle pid with respect to pidg e g calculate the potential of mean force of inserting one particle pid through an interface with center of mass pida Variant 5 prints the current value of the reaction coordinate Variant 6 prints a list of the binned values of the reaction coordinate e g dhins values between dmin and dmax Variant 7 prints the current potential of mean force for all values of the reaction coordinate considered Variant 8 prints the current force norm rather than vector for all values of the reaction coordinate considered Variant 9 loads a previous metadynamics sampling by reading a Tcl list of the potential of mean force and applied force This is especially useful to restart a simulation Note that the metadynamics scheme works seamlessly with the VIRTUAL_SITES feature allowing to define centers of mass of groups of particles as end points of the reaction coordinate One can therefore measure the potential of mean force of the distance between a particle and a molecule or interface The metadynamics scheme has as of now only been implemented for one processor MPI usage is not supported However one can speed up sampling by communicating the profile and force between independent simulations denoted walkers The print_ stat and load_s
55. on hk oneal ene Miata We de oe a Geen bl eco fale fe e poa reused A O E The MMM family of algorithms ena o ad a a A IS ee is Bs A te ek ae A ce eg E ace o Ge ee ne ek ss A E4 ELC E 5 Errors F Bibliography 210 210 212 214 215 216 218 223 1 Introduction ESPResSo is a simulation package designed to perform Molecular Dynamics MD and Monte Carlo MC simulations It is meant to be a universal tool for simulations of a variety of soft matter systems ESPResSo features a broad range of interaction poten tials which opens up possibilities for performing simulations using models with different levels of coarse graining It also includes modern and efficient algorithms for treatment of electrostatics P3M MMM type algorithms Maggs algorithm hydrodynamic interactions DPD Lattice Boltzmann and magnetic interactions It is designed to exploit the capabilities of parallel computational environments The program is being continuously extended to keep the pace with current developments both in the algorithms and software The kernel of ESPResSo is written in C with computational efficiency in mind In teraction between the user and the simulation engine is provided via a Tcl scripting interface This enables setup of arbitrarily complex systems which users might want to simulate in future as well as modifying simulation parameters during runtime 1 1 Guiding principles ESPResSo is a tool for performing compute
56. outperforms P3M easily Especially for systems with more than 1000 charges Of course if the system is not dense enough one will have to set the lattice spacing in a way to avoid several particles in one cell and the mesh will be very fine for not so many charges Also if you have lots of charges but your simulation should only run for a short time the initialization scheme takes too long in comparison But if you have dense systems with more than 1000 charges or simulations that run for many timesteps this method is definitely an option 209 E The MMM family of algorithms E 1 Introduction In the MMM family of algorithms for the electrostatic interaction a convergence factor approach to tackle the conditionally convergent Coulomb sum is used even the authors of the original MMM method have no idea what this acronym stands for Instead of defining the summation order one multiplies each summand by a continuous factor c B Tij kim such that the sum is absolutely convergent for 8 gt 0 but c 0 1 The energy is then defined as the limit 6 gt 0 of the sum i e PB is an artificial convergence parameter For a convergence factor of e 8 iim the limit is the same as the spherical limit and one can derive the classical Ewald method quite conveniently through this approach 48 To derive the formulas for MMM one has to use a different convergence factor namely e 4l ii em which defines the alternative energy N 7 li
57. points per side with which to divide the bilayer for height profile analyses e str 4 0 Distance of a tail bead from bilayer midplane beyond which a lipid is deemed to be a stray lipid Sets up the analysis package for a simulation run or analysis run that is about to be performed This routine needs to be called before any analysis can be performed mbtools analysis do_analysis 162 Calls all of the analyze routines corresponding to commands setup in setup_analysis do_analysis should be called only after setup_analysis has already been called mbtools analysis reset_averages Calls all of the resetav routines corresponding to commands setup in setup_analysis These routines vary from command to command but they typically reset the storage and counter variables used for analysis results reset_averages is typically only called internally by print_averages mbtools analysis print_averages Calls all of the printav routines corresponding to commands setup in setup_analysis These routines typically print results to a file buffer After printing the reset_averages routine is called internally print_averages should be called only after setup_analysis has already been called 14 4 2 Available analysis routines boxl verbose output time_vs_boxl Simply obtains the box dimensions from ESPResSo clusters alipid arg verbose output time_vs_clust sizehisto format 05d time e alipid 1 29 Value
58. potential is 0 This potential allows capping the force using inter forcecap see section 5 9 5 5 1 9 Soft sphere interaction Syntax inter typel type soft sphere a Nn Teut Toffset Required features SOFT_SPHERE Description This defines a soft sphere interaction between particles of the types typel and type2 which is defined by a single power law V r a r ae 5 9 for r lt Teut and V r 0 above There is no shift implemented currently which means that the potential is discontinuous at r Teut Therefore energy calculations should be used with great caution 5 1 10 Hat interaction Syntax inter typel type hat Fmax Ye Required features HAT Description This defines a simple force ramp between particles of the types type1 and type2 The maximal force Fmax acts at zero distance and zero force is applied at distances re and bigger For distances smaller than re the force is given by F r Fnax 1 z 5 10 Te for distances exceeding re the force is zero The potential energy is given by Pret r EE 1 5 11 2re which is zero for distances bigger than re and continuous at distance re This is the standard conservative DPD potential and can be used in combination with inter DPD 5 9 2 The potential is also useful for live demonstrations where a big time step may be employed to obtain quick results on a weak machine for which the physics do not need to be entirely correct 48 5 1 1
59. quantities Syntax 1 correlation id write_to_file filename 2 correlation id print 3a correlation id print averagel variance1 correlation_time 3b correlation id print average_errorbars Description Variant 1 writes the current status of the correlation estimate to the specified filename If the file exists its contents will be overwritten Output format The output looks as follows taul n_samples C1 C2 Cn tau2 n_samples C1 C2 Cn 123 Where each line corresponds to a given value of tau n_samples is the number of samples which contributed to the correlation at this level and C are the individual components of the correlation Variant 2 returns the current status of the correlation estimate as a Tcl variable Output format The output looks as follows taul n_samples C1 C2 Cn tau2 n_samples Ci C2 Cn Variants 3a and 3b return the corresponding estimate of the statistical property as a Tcl variable averagel prints the average of observablel variancel prints the variance of observablel correlation_time prints the estimate of the correlation time average_errorbars prints the estimate of the error of the average based on the method according to same as used by the uwerr command 9 2 6 The correlation algorithm multiple tau correlator Here we briefly describe the multiple tau correlator which is implemented in ESPResSo For a more detailed description and discussion of its behav
60. r_cut w and w are identical to 0 in both cases The friction dissipative and noise random term are coupled via the fluctuation dissipation theorem The friction term is a function of the relative velocity of particle pairs The DPD thermostat is better for dynamics than the Langevin thermostat since it mimics hydrodynamics in the system i When using a Lennard Jones interaction r_cut 260 is a good value to choose so that the thermostat acts on the relative velocities between nearest neighbor particles Larger cutoffs including next nearest neighbors or even more are unphysical gamma is basically an inverse timescale on which the system thermally equilibrates Values between 0 1 and 1 are o k but you propably want to try this out yourself to get a feeling for how fast temperature jumps during a simulation are The dpd thermostat does not act on the system center of mass motion Therefore before using dpd you have to stop the center of mass motion of your system which you can achieve by using the command galilei_transform 6 9 This may be repeated once in a while for long runs due to round off errors check this with the command system_CMS_velocity Two restrictions apply for the dpd implementation of ESPResSo e As soon as at least one of the two interacting particles is fixed see 4 on how to fix a particle in space the dissipative and the stochastic force part is set to zero for both particles you should only change this hardcod
61. rdf lt rdf gt part_type_list_a part_type_list_b rmin rmax rbins analyze structurefactor type order analyze vanhove type rmin rmax rbins tmax analyze centermass partype analyze momentofinertiamatrix typeid gpb Manningparameter outercellradius innercellradius Riiie I SISIERELSL Ee fe ofo Rio analyze find_principal_axis typeid analyze gyration_tensor typeid analyze aggregation dist_criteria s_mol_id f_mol_id ER o Re o E min_contact charge_criteria 192 analyze necklace pearl_threshold back_dist space_dist first length 102 analyze holes typeidprobe Mesh_size 102 analyze energy 102 analyze energy total kinetic coulomb magnetic analyze energy bonded bondid analyze energy nonbonded typeid1 typeid2 analyze pressure 103 analyze pressure total analyze pressure totals ideal coulomb tot_nonbonded_inter tot_nonbonded_intra analyze pressure bonded bondid analyze pressure nonbonded typeid1 typeid2 analyze pressure nonbonded_intra typeid analyze pressure nonbonded_inter typeid analyze stress_tensor 104 analyze stress_tensor total analyze stress_tensor totals ideal coulomb tot_nonbonded_inter tot_nonbonded_intra analyze stress_tensor bonded bond ype analyze stress_tensor nonbonded typeid1 typeid2 analyze stress_tensor nonbonded_intra typeid analyze stress_tensor nonbonded_inter typeid analyze local_stress_tensor periodic_x periodic_y periodic_z range_start_x
62. required to be able to compile and use ESPResSo Tcl Tk ESPResSo requires the Toolkit Command Language Tcl Tk in the version 8 3 or later Some example scripts will only work with Tcl 8 4 You do not only need the interpreter but also the header files and libraries Depending on the operating system these may come in separate development packages If you want to use a graphical user interface GUI for your simulation scripts you will also need Tk FFTW For some algorithms e g PM ESPResSo needs the FFTW library version 3 or later P for Fourier transforms Again the header files are required MPI Finally if you want to use ESPResSo in parallel you need a working MPI environ ment that implements the MPI standard version 1 2 1 6 Syntax description Throughout the user s guide formal definitions of the syntax of several Tcl commands can be found The following conventions are used in these descriptions e Different variants of a command are labeled 1 2 e Keywords and literals of the command that have to be typed exactly as given are written in typewriter font e If the command has variable arguments they are set in italicfont The descrip tion following the syntax definition should contain a detailed explanation of the argument and its type e alt1 alt2 specifies that one of the alternatives alt1 or alt2 can be used http www tcl tk http www fftw org 10 e argument specifies
63. roughly like this proc init id temp 92 create output files for temperature temp set f open out temp dat w close f init_particle_positions thermostat langevin temp 1 0 equilibration_integration global config set config id part setmd time The last two lines are only necessary if each instance of ESPResSo handles more than one configuration e g if you have 300 temperatures but only 10 ESPResSo processes i e load 30 In this case all user provided routines need to save and restore the configurations Saving the time is not necessary because the simulation tine across swaps is not meaningful anyways it is however convenient for investigating the temperature history of individual configurations A typical propagation routine accordingly looks like this proc perform id temp global config particle delete foreach p lindex config id 0 eval part p setmd time lindex config id 1 thermostat langevin temp 1 0 set f open out temp dat a integrate 1000 puts f setmd time analyze energy close f set config id part setmd time Again the saving and storing of the current particle properties in the config array are only necessary if there is more than one configuration per process In practice one will rescale the velocities at the beginning of perform to match the current temperature otherwise the thermostat needs a short time to equilibrate The energies necess
64. setmd variable value Description Variant 1 returns the value of the ESPResSo global variable variable variant 2 can be used to set the variable variable to value The in variant 2 means that for some variables more than one value can be given example setmd boxl 5 5 5 The following global variables can be set box_1 double 3 Simulation box length Note that if you change the box length during the simulation the folded particle coordinates will remain the same i e the particle stay in the same image box but at the same relative position in their image box If you want to scale the positions use the change_volume command cell_grid int 3 read only Dimension of the inner cell grid cell_size double 3 read only Box length of a cell dpd_gamma double read only Friction constant for the DPD thermostat dpd_r_cut double read only Cutoff for DPD thermostat gamma double read only Friction constant for the Langevin thermostat integ_switch int read only Internal switch which integrator to use local_box_1 int 3 read only Local simulation box length of the nodes max_cut double read only Maximal cutoff of real space interactions max_cut_nonbonded double read only Maximal cutoff of nonbonded real space interactions max_cut_bonded double read only Maximal cutoff of bonded real space interac tions max_num_cells int Maximal number of cells for the link cell algorithm Reason able v
65. several miscellaneous utility procedures for performing tasks such as warmup setting tabulated interactions designating molecules to be trapped and a variety of topology related sorting or data analysis functions The analysis part of the mbtools package is designed to wrap together all the analysis for a simulation into a single simple interface At the beginning of the simulation the user specifies which analyses should be performed by appending its name and arguments to a variable analysis_flags After the analysis is setup one can then simply call do_ analysis to perform all the specified proceedures Analysis will store a data value each time do_analysis is called Then when a call to print_averages is made the average of all stored values is printed to a file and the store of values is reset to nil This documentation was written by Ira R Cooke and published on his website It has been transcripted by Tristan Bereau 158 14 2 Installing and getting started Since mbtools is provided as part of the espresso molecular dynamics simulation package you will need to download and install Espresso before you can use it Espresso can be downloaded free from http espressomd org Once you have installed espresso you can find mbtools by looking inside the packages subdirectory Inside the packages mbtools directory you will see a directory for each of the mbtools subpackages as well as an examples directory All of the examples scripts
66. should work out of the box except those involving colloids which require you to install icover sh see documentation for hollowsphere molecule type To run the simplebilayer example cd to the examples directory and then type ESPRESSO_SOURCE PLATFORM Espresso scripts main tcl simplebilayer tcl The first part of this command is simply the full path to the appropriate espresso executable on your machine when running on multiple processors Obviously you will need to have the ESPRESSO_SOURCE and PLATFORM environment variables set for it to work After this executable the relative paths to two tcl scripts are given The first of these main tcl is given as an argument to espresso and is therefore interpreted first by the espresso tcl interpreter The second tcl script simplebilayer tcl is in turn passed as an argument to main tcl Why separate the tcl commands into two files This is really a matter of preference but if we keep all of the key commands and complex coding in a single file main tcl and delegate simple parameter setting to a separate file it tends to be much easier to manage large numbers of jobs with regularly changing requirements Regardless of your personal preferences the important point to note is that all of the important commands are contained in main tcl and you should probably start there to get an understanding for how mbtools works Running the simplebilayer example should produce a directory called simplebilayer which co
67. such system specification in turn should be a list consisting of a geometry and a list detailing the number of each molecule type i e set system_spec geometry n_molslist The geometry should be specified as a list with two elements The first element should be a string geometry identifying this list as a geometry The second element is a string containing the name of a geometry type mygeometry followed by arguments to be passed to the routine create_mygeometry The n_molslist should be specified as a list with two elements The first element should be a string n_molslist identifying this list as an n_molslist The second element is a list each element of which specifies a molecule type and the number of such molecules e bozxl A list containing the lengths of each of the box side lengths e moltypes A list each element of which specifies a molecule type and type informa tion The exact format and requirements of this list are detailed for each molecule separately see below for a list of molecule types and their requirements however regardless of mol type the first two elements of the list must be a moltypeid and a string specifying the moltype respectively Sets up the system including generating topologies and placing molecules into specified geometries Each geometry and list of molecules to be placed into that geometry are grouped into a system spec Example The following code sets out the molecule types to be used
68. the Langevin thermostat as documented in section 6 2 sets the temperature T individually for the particle with id pid This allows to simulate systems containing particles of different temperatures Caution this has no influence on any other thermostat then the Langevin thermostat e gamma g If used in combination with the Langevin thermostat as documented in section 6 2 sets the frictional coefficient T individually for the particle with id pid This allows to simulate systems containing particles with different diffusion constants Caution this has no influence on any other thermostat then the Langevin thermostat e rotation rot Specifies whether a particle s rotational degrees of freedom are integrated value of 1 or not 0 If set to zero the content of the torque and omega variables are meaningless The default is 1 4 1 2 Getting particle properties Syntax 1 part pid print id pos type folded_position type q v f fix ext_force ext_torque bond connections range 2 part Description Variant 1 will return a list of the specified properties of particle pid or all properties if no keyword is specified Variant 2 will return a list of all properties of all particles Example part 40 print id pos q bonds will return a list like 40 8 849 1 8172 1 4677 1 0 This routine is primarily intended for effective use in Tcl scripts When the keyword connection is specified it returns t
69. the outward pointing normal vectors for every surface element The parameter convergence_criterion allows to specify the accuracy of the iteration It corresponds to the maximum relative change of any of the interface particle s charge After max_iterations the iteration stops anyways The dielectric constant in bulk i e outside the dielectric walls is specified by eps_out A homogenous electric field can be added to the calculation of dielectric boundary forces by specifying it in the parameter ext_field 67 Quick setup of dielectric interfaces Syntax 1 dielectric sphere center cx cy cz radius r res res 2 dielectric wall normal nz ny nz dist d res res 3 dielectric cylinder center cx cy cz axis av ay az radius r direction d 4 dielectric pore center cx cy cz axis ax ay az radius r length l smoothing radius rs res res Description The command dielectric allows to conveniently create dielectric interfaces similar to the constraint and the Ibboundary command Currently the creation of spherical cylin drical and planar geometries as well as a pore geometry is supported It is implemented in Tcl and places particles in the right positions and adds the correct values to the global Tcl variables icc_areas icc_normals icc_sigmas icc_epsilons and increases the global Tcl variable varn_induced_charges Thus after setting up the shapes it is still necessary to register them by calling iccp3m usually in the following way iccp3m n_in
70. to read in the information using ESPResSo note that interactions must be stored before particles before bonding information as for the bonds to be set all particles and all interactions must already be known to ESPResSo set in open gzip cd checkpoint block gz r while blockfile in read auto eof close in This is basically all you need to restore the information in the blockfile overwriting the current settings in ESPResSo 131 10 3 Writing VTF files The formats VTF WTF Trajectory Format VSF VTF Structure Format and VCF VTF Coordinate Format are formats for the visualization software VMDBof They are intended to be human readable and easy to produce automatically and modify The format distinguishes between structure blocks that contain the topological infor mation of the system i e the system size particle names types radii and bonding information amongst others while coordinate blocks a k a as timestep blocks contain the coordinates for the particles at a single timestep For a visualization with VMD one structure block and at least one coordinate block is required Files in the VSF format contain a single structure block files in the VCF format contain at least one coordinate block while files in the VTF format contain a single structure block first and an arbitrary number of coordinate blocks afterwards thus allowing to store all information for a whole simulation in a single file For
71. y z Sets the rotational velocity of this particle in the particle s co rotating frame torque x y z Sets the rotational force of this particle in global frame When printing the values using part id_particle print there is an alternative named tbf which gives you the values of the torque in the body frame Be aware the values obtained when printing using torque are computed in the frame laboratory and are the ones one should usually look at Nonetheless in case you introduce torques using torque option espresso will assume they are given in the body frame Thus tbf is useful to know which should be the numerical values you should reintroduce in order to have exactly the same conformation Again any torque set this way is lost during the next integration step e rinertia z y 2 Sets the diagonal elements of this particles rotational inertia tensor These correspond with the inertial moments along the coordinate axes in the particle s co rotating coordinate system When the particle s quaternions are set to 1 0 0 0 the co rotating and the fixed lab frame are co aligned fix x y z Fixes the particle in space By supplying a set of 3 integers as ar guments it is possible to fix motion in x y or z coordinates independently For example fiz 0 0 1 will fix motion only in z Note that fiz without arguments is equivalent to fiz 1 1 1 unfix Release any external influence from the particle ext_force xz y z An additi
72. 0 Lennard Jones interaction with angular dependence between particles of the types type1 and type2 These two particles need two bonded partners oriented in a symmetric way They define an orientation for the central particle The purpose of using bonded partners is to avoid dealing with torques therefore the interaction does not need the ROTATION feature The angular part of the potential minimizes the system when the two central beads are oriented along the vector formed by these two particles The shaded beads on the image are virtual particles that are formed from the orientation of the bonded partners connected to the central beads They are used to define angles The potential is of the form 12 10 U Tik Ojik Oikn 5 6 cos Ojik cos likn 5 14 where rj is the distance between the two central beads and each angle defines the orientation between the direction of a central bead determined from the two bonded partners and the vector rik Note that the potential is turned off if one of the angle is more than 7 2 This way we don t end up creating a minimum for an anti parallel configuration Unfortunately the bonded partners are not seeked dynamically One has to keep track of the relative positions of the particle IDs This can be done by setting the parameters bla blb 624 and b2 Say the first bead type1 has particle ID n then one should set the simulation such as its two bonded partners have particle IDs n bl
73. 0 to the energy of a charge qi at position xi Yi zi in the central layer Since z z lt nAz we have z zi nAz NA zj zi and z zi n z nA Zj and hence the interaction energy from those two image layers with the charge q vanishes by charge neutrality N N Zig Y 05 l2 zi HnAgl z zi nAg 4rqind Y 0 0 j 1 j l 215 The only errors occurring are those coming from the approximation of assuming ho mogeneously charged infinite sheets instead of discrete charges This assumption should become better when increasing the distance nA from the central layer However in a naive implementation even large gap sizes will result in large errors This is due to the order of summation for the standard Ewald sum which is spherical while the above approach orders the cells in layers called slab wise summation Smith has shown that by adding to the Ewald energy the term 27M E 21M where M is the total dipole moment one obtains the result of a slab wise summation instead of the spherical limit 48 Although this is a major change in the summation order the difference is a very simple term In fact Smith shows that changes of the summation order always result in a difference that depends only on the total dipole moment Using the far formula of MMM2D one can calculate the contributions of the addi tional layers up to arbitrarily precision even for small gap s
74. 0000 particles and a typical cell grid size of 20 these two cell layers consume roughly 200 KBytes which nearly fits into the L2 cache Therefore every cell has to be read from the main memory only once per force calculation 183 16 Getting involved Up to date information about the development of ESPResSo can be found at the web page http espressomd org As the important information can change in time we will not describe its contents in detail but rather request the reader to go directly to the URL Among other things one can find information about the following topics there e FAQ e Latest stable release of ESPResSo and older releases Obtaining development version of ESPResSo e Archives of both developers and users mailing lists e Registering to ESPResSo mailing lists Submitting a bug report 16 1 Community support and mailing lists If you have any questions concerning ESPResSo which you cannot resolve by yourself you may post a message to the mailing list Instructions on how to register to the mailing lists and post messages can be found on the homepage http espressomd org Before posting a question and waiting for someone to answer it may be useful to search the mailing list archives or FAQ and see if you can get the answer immediately For several reasons it is recommended to send all questions to the mailing lists rather than to contact individual developers e All registered users get your messa
75. 1 Hertzian interaction Syntax inter typel type2 hertzian o Required features HERTZIAN Description This defines an interaction according to the Hertzian potential between particles of the types typel and type2 The Hertzian potential is defined by NE ge ie vin fel 5 es 5 12 The potential has no singularity and is defined everywhere the potential has nondiffer entiable maximum at r 0 where the force is undefined 5 1 12 Gaussian Syntax inter typel type gaussian O Teut Required features GAUSSIAN Description This defines an interaction according to the Gaussian potential between particles of the typers typel and type2 The Gaussian potential is defined by V r aed r lt fout 5 13 0 T gt Teut The Gaussian potential is smooth except at the cutoff and has a finite overlap energy of e It can be used to model e g overlapping polymer coils Currently there is no shift implemented which means that the potential is discontin uous at r reyt Therefore use caution when performing energy calculations However you can often choose the cutoff such that the energy difference at the cutoff is less than a desired accuracy since the potential decays very rapidly 5 2 Anisotropic non bonded interactions 5 2 1 Directional Lennard Jones interaction Syntax inter typel type lj angle O Tout bla bl b2q b2 fcap 20 02 hi Required features LJ_ANGLE 49 Description Specifies a 12 1
76. 1 returns all stored configurations while variant 2 returns only the number of stored configurations 112 Output format variant 1 x1 yl 21 2 92 22 a 8 4 uwerr Computing statistical errors in time series Syntax 1 uwerr data nrep col s_tau plot 2 uwerr data nrep f s_tau f_args plot Description Calculates the mean value the error and the error of the error for an arbitrary numerical time series according to Wolff 57 Arguments e data is a matrix filled with the primary estimates az from R replica with N1 Na Nr measurements each 1 1 1 1 1 1 1 ag 3 2 1 2 1 2 1 a a5 a3 data ae a ait 1 2 1 2 1 2 ay a2 3 Nr R Np R _Np R a R Ay R a3 R e nrep is a vector whose elements specify the length of the individual replica nrep Ni N2 NR ef is a user defined Tcl function returning a double with first argument a vector which has as many entries as data has columns If f is given instead of the column the corresponding derived quantity is analyzed ef_args are further arguments to f es tau is the estimate S 7 Tin as explained in section 3 3 of 57 The default is 1 5 and it is never taken larger than min N 2 e plot If plot is specified you will get the plots of T T 0 and Tint vs W The data and gnuplot script is written to the current directory 113 Output format mean error error_of_error act error_of_act Q where act denotes the integrated autocorr
77. 1 0 cnt this is necessary since 1 cnt is interpreted 17 g r N Figure 2 2 Radial distribution functions gy r between equal charges rectangles and g r for opposite charges circles The plus symbols denote g r for an uncharged system as an integer division which results in 0 for cnt gt 1 argv is a predefined variable it contains all the command line parameters Therefore this script should be called like Espresso script config The printing of the calculated radial distribution functions is simple Add to the end of the previous snippet the following lines set plot open rdf data w puts plot r rdf r foreach r rlist rdf avg_rdf 4 puts plot r rdf close plot This instructs the Tcl interpreter to write the avg_rdf to the file rdf data in gnuplot compatible format Fig shows the resulting radial distribution functions averaged over 100 configurations In addition the distribution for a neutral system is given which can be obtained from our simulation script by simply removing the command inter coulomb and therefore not turning on P M The code example given before is still quite simple and the reader is encouraged to try to extend the example a little bit e g by using differently sized particle or changing the interactions If something does not work ESPResSo will give comprehensive error messages which should make it easy to identify mistakes
78. 1 1 Source and build directories Usually when a program is compiled the resulting binary files are put into the same directory as the sources of the program In ESPResSo s build system the source directory that contains all the source files is completely separated from the build directory where the files created by the build process are put The location of the build directory is the current working directory at the time when configure is called In this way you can build several variants of ESPResSo each variant having different activated features and for as many platforms as you want All further commands concerning compiling and running ESPResSo have to be called from the build directory Example When the source directory is srcdir i e the files where unpacked to this directory then the build directory can be set to builddir by calling the configure script from there cd builddir srcdir configure make Espresso 3 1 2 Options The behaviour of configure can be controlled by the means of command line options In the following only those command line options that are specific to ESPResSo will be explained For a complete list of options and explanations thereof call configure help with myconfig MYCONFIG_HEADER This option sets the name of the local configura tion header see 3 4 on page 24 It defaults to myconfig h with mpi yes no guess without mpi By default configure will auto matica
79. 2 2 W a y z Io exp i w w w Equation is a generalization of the formula presented by H fling et al 24 For more information see references therein Per each 3 dimen sions of the observable one dimension of the correlation output is produced If fcs_acf is used with other observables than particle_positions the physical meaning of the result is unclear 121 dt The time interval of sampling data points When autoupdate is used dt has to be a multiple of timestep It is also used to produce time axis in real units Warning if dt is close to the timestep autoupdate is strongly recommended Otherwise cpu time is wasted on passing the control between the script and kernel tau_max This is the maximum value of 7 for which the correlation should be computed Warning Unless you are using the multiple tau correlator choosing tau_maz of more than 100dt will result in a huge computational overhead In a multiple tau correlator with reasonable parameters tau_max can span the entire simulation without too much additional cpu time tau_lin The number of data points for which the results are linearly spaced in tau This is a parameter of the multiple tau correlator If you want to use it make sure that you know how it works By default it is set equal to tau mags which results in the trivial linear correlator By setting tau_lin tau_max the multiple tau correlator is switched on In many cases tau_lin 16 is a good choi
80. 2d slab geometries Comput Phys Commun 148 3 327 348 2002 5 7 3 Axel Arnold Jason de Joannis and Christian Holm Electrostatics in Periodic Slab Geometries I II J Chem Phys 117 2496 2512 2002 5 7 6 Axel Arnold Olaf Lenz Stefan Kesselheim Rudolf Weeber Florian Fahrenberger Dominic Roehm Peter Kosovan and Christian Holm ESPResSo 3 1 molecu lar dynamics software for coarse grained models In Michael Griebel Christian Rieger and Marc Alexander Schweitzer editors Proceedings of the Sixth Interna tional Workshop on Meshfree Methods for Partial Differential Equations Lecture Notes in Computational Science and Engineering Springer Berlin Germany sub mitted V Ballenegger A Arnold and J J Cerda Simulations of non neutral slab sys tems with long range electrostatic interactions in two dimensional periodic bound ary conditions J Chem Phys 131 9 094107 2009 doi 10 1063 1 3216473 URL http link aip org link JCP 131 094107 1 A Brodka Ewald summation method with electrostatic layer correction for in teractions of point dipoles in slab geometry Chem Phys Lett 400 62 67 2004 9 8 2 Juan J Cerda V Ballenegger O Lenz and C Holm P3M algorithm for dipolar interactions J Chem Phys 129 234104 2008 5 7 1 5 8 1 218 11 12 13 14 15 16 17 18 19 20 21 22 23 24 I Cimr k M Gusenbauer and T Schrefl Mod
81. 9 LANGEVIN PER_PARTICLE 1 ROTATION_PER_PARTICLE je part pid print id pos type folded_position type q v 28 fix ext_force ext_torque bond connections range part part pid delete 29 part deleteall part auto_exclusions range 29 part delete_exclusions Required features EXCLUSIONS polymer num_polymers monomers_per_chain bond_length 30 start pid pos x y z mode RW SAW PSAW shield trymax charge q distance denargoa types typeidneutral tYPCidenarged bond bondid angle 0 x y 2 constraints Required features ELECTROSTATICS 2CONSTRAINTS counterions N start pid mode SAW RW shield trYmax 32 charge val type typeid Required features ELECTROSTATICS salt N N_ start pid mode SAW RW shield trymax 32 charges val val_ t types typeid typeid_ rad r Required features ELECTROSTATICS diamond a bond_length monomers_per_chain counterions Nc 33 charges Valnode Valmonomer vale distance dal nonet Required features ELECTROSTATICS 186 icosaeder a monomers_per_chain counterions Nc 34 1 charges Valmonomers Valcr distance dcharged Required features ELECTROSTATICS crosslink num_polymer monomers_per_chain start pid catch eaten 35 distLink link_dist distChain chain_dist FENE bondid trials trymax copy_particles set idi id2 range from to shift sx s
82. CTROSTATICS dielectric sphere center cx cy cz radius r res res dielectric wall normal nz ny nz dist d res res dielectric cylinder center cx cy cz axis ax ay az radius r direction d dielectric pore center cx cy cz axis ax ay az radius r length l smoothing_radius rs res res J 3 13 oo ep el s S z z z Aa z gt E inter magnetic 0 0 67 inter magnetic inter magnetic parameters inter magnetic lg p3m Trax mesh cao alpha Required features DIPOLES inter magnetic lg p3m tune tunev2 accuracy accuracy r_cut Test mesh mesh cao cao alpha a Required features DIPOLES inter magnetic mdlc accuracy gap_size far_cutoff Required features DIPOLES inter magnetic lg dawaanr Required features DIPOLES inter magnetic lg mdds n_cut value_n_cut Required features DIPOLES MAGNETIC_DIPOLAR_DIRECT_SUM inter typel type2 tunable_slip T YL Teut t Ur Vy vz 70 Required features TUNABLE_SLIP inter typel type2 inter_dpd gamma r_cut wf tgamma tr_cut twf 70 Required features INTER_DPD inter typeid1 typeid1 comfixed flag 71 Required features COMFIXED inter typeid1 typeid2 comforce flag dir force fratio 72 Required features COMFORCE inter forcecap Fmax individual 72 setmd variable setmd variable value thermostat 73 thermostat off thermostat parameters thermostat langevin temperature gamma thermostat ghmc temperature n
83. D simulation here The code pieces can be copied step by step into a file which then can be run using Espresso file from the ESPResSo source directory Our script starts with setting up the initial configuration Most conveniently one would like to specify the density and the number of particles of the system as parameters set n_part 200 set density 0 7 set box_1 expr pow n_part density 1 3 These variables do not change anything in the simulation engine but are just standard Tcl variables they are used to increase the readability and flexibility of the script The box length is not a parameter of this simulation it is calculated from the number of particles and the system density This allows to change the parameters later easily e g to simulate a bigger system The parameters of the simulation engine are modified by the setmd command For example http www tcl tk man tcl8 5 tutorial tcltutorial html 13 setmd box_1 box_1 box_1 box_1 setmd periodic 1 1 1 defines a cubic simulation box of size box_1 and periodic boundary conditions in all spatial dimensions We now fill this simulation box with particles set q 1 set type 0 for set i 0 i lt n_part incr i set posx expr box_1 t_random set posy expr box_1 t_random set posz expr box_1 t_random set q expr q set type expr 1 type part i pos posx posy posz q q type type This loop adds n_part particles at random positi
84. E ROTATION_PER_PARTICLE Description This command modifies particle data namely position type monomer ion charge velocity force and bonds Multiple properties can be changed at once If you add a new particle the position has to be set first because of the spatial decomposition Arguments e pid elpos x y z Sets the position of this particle to x y z e type typeid Restrictions typeid gt 0 The typeid is used in the inter command see section 5 on page 43 to define the parameters of the non bonded interactions between different kinds of particles elv vx vy vz Sets the velocity of this particle to vx vy vz The velocity remains variable and will be changed during integration elf fx fy fz Set the force acting on this particle to fx fy fz The force remains variable and will be changed during integration However whereas the velocity is modified with respect to the velocity you set upon integration the force it recomputed during the integration step and any force set in this way is lost during the integration step e bond bondid pid2 Restrictions bondid gt 0 pid2 must be an existing parti cle The bondid is used for the inter command to define bonded interactions 26 ebond delete Will delete all bonds attached to this particle q charge Sets the charge of this particle to q quat ql q2 q3 q4 Sets the quaternion representation of the rotational position of this particle omega x
85. EXCLUSIONS Description Variant 1 will create exclusions for all particles pairs connected by not more than range bonds range defaults to 2 This is typically used in atomistic simulations where nearest and next nearest neighbour interactions along the chain have to be omitted since 29 they are included in the bonding potentials For example if the system contains particles 0 100 where particle n is bonded to particle n 1 for 1 lt n lt 100 then it will result in the exclusions e particle 1 does not interact with particles 2 and 3 e particle 2 does not interact with particles 1 3 and 4 e particle 3 does not interact with particles 1 2 4 and 5 e Variant 2 deletes all exclusions currently present in the system 4 2 Creating groups of particle 4 2 1 polymer Setting up polymer chains Syntax polymer num_polymers monomers_per_chain bond_length start pid pos x y z mode RW SAW PSAW shield trymax charge q distance del types typeidneutral typetdcharged bond bondid angle 9 x y 2 constraints Required features ELECTROSTATICS CONSTRAINTS Description This command will create num_polymers polymer or polyelectrolyte chains with monomers _per_chain monomers per chain The length of the bond between two adjacent monomers will be set up to be bond_length Arguments enum_polymers Sets the number of polymer chains e monomers_per_chain Sets the number of monomers per
86. MASS_LIN DPD_MASS_RED ELECTROSTATICS ESK DEBUG ESR DEBUG EVENT_DEBUG EWALD_DEBUG EXCLUSIONS EXTERNAL_FORCES FENE_DEBUG FFT_DEBUG FORCE_CORE FORCE_DEBUG GAY BERNE GHOST_DEBUG GHOST_FORCE_DEBUG GRID_DEBUG HALO_DEBUG HERTZIAN INTEG_DEBUG INTER_DPD 199 INTER_RF LANGEVIN PER PARTICLE LANGEVIN PER PARTICLE LATTICE DEBUG LB 199 LB DEBUG LB_ELECTROHYDRODYNAMICS 199 LENNARD JONES LENNARD_JONES_GENERIC 200 LJ_ANGLE LJ_DEBUG LJ_WARN_WHEN_CLOSE LJCOS LJCOS2 MAGGS_DEBUG MASS MEM_DEBUG METADYNAMICS MODES MOL CUT MOLFORCES MOLFORCES_DEBUG MORSE MORSE_DEBUG MPLCORE NEMD NO_INTRA_NB 200 NPT 199 OLD_RW_VERSION 224 OLD_DIHEDRAL ONEPART DEBUG OVERLAPPED 199 P3M_DEBUG PARTIAL PERIODIC PARTICLE_DEBUG POLY DEBUG PTENSOR DEBUG RANDOM DEBUG ROTATION ROTATION PER PARTICLE ROTATIONAL INERTIA SMOOTH STEP SOFT_SPHERE STAT_DEBUG TABULATED 199 THERMO_DEBUG THERMOSTAT IGNORE_NON VIR TUAL TRANS_DPD 199 TUNABLE SLIP VERLET DEBUG VIRTUAL SITES COM VIRTUAL _SITES_DEBUG VIRTUAL SITES NO_VELOCITY 198 VIRTUAL SITES _RELATIVE 198 VIRTUAL SITES THERMOSTAT 198 FENE bond FFTW finding holes form factor of a chain 110 fsi Tcl command gamma global variable Gaussian interaction Gay Berne interaction Generic Lennard Jones interaction g
87. NGEVIN_PER_PARTICLE Allows to define the temperature and friction coefficient for individual particles See for details 198 e CATALYTIC_REACTIONS Allows the user to define three particle types to be reactant catalyzer and product Reactants get converted into products in the vicinity of a catalyst according to a used defined reaction rate constant It is also possible to set up a chemical equilibrium reaction between the reactants and products with another rate constant See section 6 8 for details e OVERLAPPED e COLLISION_DETECTION Allows particles to be bound oo collision See section e OLD_RW_VERSION This switches back to the old wrong random walk code of the polymer Only use this if you rely on the old behaviour and know what you are doing In addition there are switches that enable additional features in the integrator or thermostat e NEMD Enables the non equilbrium shear MD support see section 6 3 on page 80 e NPT Enables an on the fly NPT integration scheme see section 6 2 4 on page 80 e DPD Enables the dissipative particle dynamics thermostat see section pase 77 e TRANS_DPD Enables the transversal dissipative particle dynamics thermostat see section 6 2 3 on page 79 e INTER_DPD Enables the dissipative particle dynamics thermostat implemented as an interaction allowing to choose different parameters between different particle types see section 6 2 3 on page 79 e INTER_RF e DPD_MASS
88. TROSTATICS Description MMM2D coulomb method for systems with periodicity 1 1 0 Needs the layered cell system The performance of the method depends on the number of slices of the cell system which has to be tuned manually It is automatically ensured that the maximal pairwise error is smaller than the given bound The far cutoff setting should only be used for testing reasons otherwise you are more safe with the automatical tuning If you even don t know what it is do not even think of touching the far cutoff For details on the MMM family of algorithms refer to appendix E on page 210 The last two mutually exclusive arguments dielectric and dielectric constants allow to specify dielectric contrasts at the upper and lower boundaries of the simulation box The first form specifies the respective dielectric constants in the media which however is only used to calculate the contrasts That is specifying m p const is always identical to amp m amp 1 The second form specifies only the dielectric contrasts at the boundaries that is A ane and A ene Using this form allows to choose A 1 corresponding to metallic boundary conditions 5 7 4 MMM1D Please cite 3 BIBTEX key mmmid in file doc ug citations bib when using MMMI1D Syntax 1 inter coulomb lg mmmid switch radius bessel_cutoff mazimal_pairwise_error 2 inter coulomb lg mmmid tune mazimal_pairwise_err
89. This defaults to 12000 e inode Specifies how many systems to run per ESPResSo instance If this is more than 1 it is the user s responsibility to manage the storage of configurations see below for examples This defaults to 1 e Rreset Specifies after how many swap trial rounds to reset the counters for the acceptance rate statistics This defaults to 10 e info specifies which output the parallel tempering code should produce 91 none parallel tempering will be totally quiet except for fatal errors comm information on client activities such as connecting is printed to stderr all print lots of information on swap energies and probabilities to stdout This is useful for debugging and quickly checking the acceptance rates This defaults to all Introduction The basic idea of parallel tempering is to run N simulations with configurations C in parallel at different temperatures T lt Ty lt lt Ty and exchange configurations between neighboring temperatures This is done according to the Boltzmann rule t e the swap probability for two configurations A and B at two different parameters T and T is given by min 1 exp 8 72 Ua Ta 8 T UB Ti B 11 Ua T B Z2 UB T2 7 1 where Uc T denotes the potential energy of configuration C at parameter T and T the corresponding inverse temperature If T is the temperature Uc is indepedent of T and T 1 kgT In this case the swap probability red
90. _RED Enables masses in DPD using reduced dimensionless mass units e DPD_MASS_LIN Enables masses in DPD using absolute mass units e LB Enables the lattice Boltzmann fluid code see section 12 on page 148 e LB_ELECTROHYDRODYNAMICS Enables the implicit calculation of electro hydrodynamics for charged particles and salt ions in an electric field B 2 Interactions The following switches turn on various short ranged interactions see section pago 23 e TABULATED Enable support for user defined interactions 199 LENNARD_JONES Enable the Lennard Jones potential LENNARD_JONES_GENERIC Enable the generic Lennard Jones potential with config urable exponents and individual prefactors for the two terms LJCOS Enable the Lennard Jones potential with a cosine tail LJCOS2 Same as LJCOS but using a slightly different way of smoothing the con nection to 0 LJ_ANGLE Enable the directional Lennard Jones potential GAY_BERNE HERTZIAN MOL_CUT NO_INTRA_NB MORSE Enable the Morse potential BUCKINGHAM Enable the Buckingham potential SOFT_SPHERE Enable the soft sphere potential SMOOTH_STEP Enable the smooth step potential a step potential with two length scales BMHTF_NACL Enable the Born Meyer Huggins Tosi Fumi potential which can be used to model salt melts Some of the short range interactions have additional features e LJ_WARN_WHEN_CLOSE This adds an additional check to the Lennard Jones po tential
91. _md phi no_flip flip random_flip 76 no_scale scale thermostat dpd temperature gamma r_cut WF wf tgamma tr_cut TWF twf 76 Required features DPD or TRANS_DPD thermostat inter_dpd temperature 77 Required features INTER_DPD thermostat npt_isotropic temperature gamma0 gammaV 77 Required features NPT 190 3 nemd exchange n_slabs n_exchange nemd shearrate n_slabs shearrate nemd off nemd nemd profile nemd viscosity Required features NEMD cellsystem domain_decomposition no_verlet_list cellsystem nsquare cellsystem layered n_layers adress set topo kind width width hybrid_width center x Rx wf wf Required features ADRESSO cuda list cuda setdevice id cuda getdevice on_collision on_collision off on_collision exception bind_centers d bondi on_collision exception bind_at_point_of_collision d bond1 bond2 type reaction reactant_type rt catalyzer_type ct product_type pt range r ct_rate k_ct eq_rate k_eq react_once on off reaction off reaction print Required features CATALYTIC_REACTIONS The current implementation also requires the use of verlet lists and domain decomposi tion kill_particle_motion rotation Required features 1 ROTATION kill_particle_forces torques Required features 1 ROTATION system_CMS system_CMS_velocity galilei_transform integrate steps integrate set nvt integrate set npt_isotropic Pert piston x y 2 cubic_box change_vol
92. _y s_z 35 constraint wall normal ng Ny n dist d type id penetrable flag 36 reflecting flag constraint sphere center Cr Cy Cz radius rad direction direction type id penetrable flag reflecting flag constraint cylinder center Cr Cy Cz axis Nz Ny Nn radius rad length length direction direction type id penetrable flag reflecting flag constraint rhomboid corner Pz Py Pz A Ax Qy dz b by by bz C Cz Cy Cz direction direction type id penetrable flag reflecting flag constraint maze nsphere n dim d sphrad r cylrad re type id penetrable flag constraint pore center Cr Cy Cz axis Ng Ny nz radius rad length length type id constraint rod center C Cy lambda lambda A constraint plate height h sigma sigma constraint ext_magn_field fr fy fz 2 3 constraint plane cell xz y z type id constraint mindist_position z y z 1 Required features CONSTRAINTS ELECTROSTATICS ROTATION DIPOLES El constraint delete num constraint force n w 00 constraint num part gc type find delete status number type Required features GRANDCANONICAL i inter inter typel type2 tabulated filename Required features TABULATED i inter typel type2 lennard jones O Teut Cshiftlauto row Teap Tmin 44 Required features LENNARD_JONES Required features LENNARD_JONES_GENERIC inter typel type lj gen O Teut Cshift To 1 2 bi ba Tcaplauto 187 inter typel type lj cos O To
93. a and n bl respectively On a linear chain for example one would typically have bla 1 and bl 1 such that the central bead and its two bonded partners have position IDs n n 1 and n 1 respectively This is surely not optimized but once the simulation is set correctly the algorithm is very fast The force can be capped using inter forcecap It might turn out to be useful in some cases to keep this capping during the whole simulation This is due to the very sharp angular dependence for small distance compared to a Two beads might come very close to each other while having unfavorable angles such that the interaction is turned off Then a change in the angle might suddenly turn on the interaction and the system will blow up the potential is so steep that one would need extremely small time steps to deal with it which is not very clever for such rare events For instance when modeling hydrogen bonds N H O C one can avoid simulating hydrogens and oxygens by using this potential This comes down to implementing a HBond potential between N and C atoms 50 The optional parameter reap is the usual cap radius The four other optional pa rameters 29 62 k describe a different interaction strength e for a subset of the simulation box The box is divided through the z plane in two different regions region 1 which creates an interaction with strength e region 2 with interaction strength e The 2nd region is def
94. a positive floating point number Make sure that you know the relevance of the P3M parameters before using P3M If you are not sure read the following references 19 23 30 141 15 16 1131 10 Tuning Coulomb P3M Syntax inter coulomb lg p3m tune tunev2 accuracy accuracy r_cut Tout mesh mesh cao cao alpha a Required features ELECTROSTATICS Description It is not easy to calculate the various parameters of the P3M method such that the method provides the desired accuracy at maximum speed To simplify this ESPResSo provides a function to automatically tune the algorithm Note that for this function to work properly your system should already contain an initial configuration of charges and the correct initial box size Also note that both provided tuning algorithms work very well on homogenous charge distributions but might not achieve the requested precision for highly inhomogenous or symmetric systems For example because of the nature of the P3M algorithm systems are problematic where most charges are placed in one plane one small region or on a regular grid The function employs the analytical expression of the error estimate for the P3M method and its real space error to obtain sets of parameters that yield the desired accuracy then it measures how long it takes to compute the coulomb interaction using these parameter sets and chooses the set with the shortest run time The function will only automatically tu
95. acy gap_size far_cutoff Required features DIPOLES Description Like ELC but applied to the case of magnetic dipoles but here the accuracy is the one you wish for computing the energy far utoff is set to a value that assuming all dipoles to be as larger as the largest of the dipoles in the system the error for the energy would be smaller thant the value given by accuracy At this moment you cannot compute the accuracy for the forces or torques nonetheless usually you will have an error for forces and torques smaller than for energies Thus the error for the energies is an upper boundary to all errors in the calculations At present the program assumes that the gap without particles is along the z direction The gap size is the length along the z direction of the volume where particles are not allowed to enter As a reference for the DLC method see 9 5 8 3 Dipolar all with all and no replicas DAWAANR Syntax inter magnetic lg dawaanr Required features DIPOLES Description This interaction calculates energies and forces between dipoles by explicitly summing over all pairs For the directions in which the system is periodic as defined by setmd periodic it applies the minimum image convention i e the interaction is effectively cut off at half a box length In periodic systems this method should only be used if it is not possible to use dipolar P3M or DLC because those methods have a far better accuracy and are much f
96. ally updating the correlation estimates The automatic updates are done within the integration loop without further user intervention The update frequency is adjusted based based on the value of dt provided when defining the correlation Variant 2 is an explicit call for an instantaneous update of the correlation estimates using the current system state It is only possible to use 2 if the correlation is not being autoupdated However it is possible to use it after autoupdate has been stopped When updating by an explicit call ESPResSo does not check if the lag time between two updates corresponds the value of dt specified when creating the correlation Variant 3 correlates all data from history which are left in the buffers Once this has been done the history is lost and no further updates are possible When a new observable value is passed to a correlation level 0 of the compression buffers of the multiple tau correlator see section for details is updated immediately Higher levels are updated only when the lower level buffers are filled and there is a need to push some values one level up When the updating is stopped a number of observable values have not reached the higher level especially when tau ax is comparable to the total simulation time and if there are many compression levels In such case variant 3 is very useful If tau_maz is much shorter it does not have a big effect 9 2 5 Printing out the correlation and related
97. also the tangential component is turned around so that a bounce back motion is performed The second variant is useful for boundaries of DPD The reflection property is only activated if an interaction is defined between a particular particle and the constraint This will usually be a lennard jones interaction with e 0 but finite interaction range Variants 7 and 8 create interactions based on electrostatic interactions The cor responding force acts in direction of the normal vector of the surface and applies to all charged particles Variant 9 does not define a surface but is based on magnetic dipolar interaction with an external magnetic field It applies to all particles with a dipol moment Variant 10 is essential for the use of tunable slip boundary interactions for mi crochannel flows like the Plane Poiseuille or Plane Couette Flow Variant 11 calculates the smallest distance to all non penetrable constraints that can be repulsive wall cylinder sphere rhomboid maze pore Negative distances mean that the position is within the area that particles should not access Helpful to find initial configurations Note that constraints are not saved to checkpoints and that they have to be reset upon restarting a simulation The resulting surface in variant 1 is a plane defined by the normal vector nz ny nz and the distance d from the origin The force acts in direction of the normal Note that the d describes the distance
98. alues are between 125 and 1000 or for some problems n otal articles Mn odes 74 max_part int read only Maximal identity of a particle This is in general not related to the number of particles max_range double read only Maximal range of real space interactions maz ut skin max_skin double read only Maximal skin to be used for the link cell verlet algo rithm This is the minimum of cell_size max_range min_global_cut double Minimal total cutoff for real space Effectively this plus the skin is the minimally possible cell size Espresso typically determines this value automatically but some algorithms e g virtual sites require you to specify it manually min_num_cells int Minimal number of cells for the link cell algorithm Reason able values range in 1076 N to 1077 N In general just make sure that the Verlet lists are not incredibly large By default the minimum is 0 but for the automatic P3M tuning it may be wise to set larger values for high particle numbers n_layers int read only Number of layers in cell structure LAYERED see sec tion 6 4 on page 81 n_nodes int read only Number of nodes n_part int read only Total number of particles n_part_types int read only Number of particle types that were used so far in the inter command see chaptertcl inter node_grid int 3 3D node grid for real space domain decomposition optional if unset an optimal set is chosen automatical
99. alysis desired afterwards Note that the time at which configurations were taken is not stored The most observ ables that work with the set of stored configurations do expect that the configurations are taken at equidistant timesteps Note also that the stored configurations can be written to a file and read from it via the blockfile command see section 10 2 on page 127 8 3 1 Storing and removing configurations analyze append analyze remove indez 3 analyze replace index analyze push size 5 analyze configs config Description Variant 1 appends the current configuration to the set of stored configurations Variant 2 removes the indexth stored configuration or all if index is not specified Variant 3 will replace the indexth configuration with the current configuration Variant 4 will append the current configuration to the set of stored configuration and remove configurations from the beginning of the set until the number of stored configurations is equal to size If size is not specified only the first configuration in the set is removed Variants 1 to 4 return the number of currently stored configurations Variant 5 will append the configuration config to the set of stored configurations config has to define coordinates for all configurations in the format el y1 z1 22 y2 22 2 8 3 2 Getting the stored configurations Syntax 1 analyze configs 2 analyze stored Description Variant
100. and irreversible adhesion processes Two methods of binding are available e bind_centers adds a bonded interaction between the colliding particles at the first collision This leads to the distance between the particles being fixed the particles can however still slide around each other The parameters are as follows d is the distance at which the bond is created bond1 denotes a pair bond and is the type of the bond created between the colliding particles Particles that are already bound by a bond of this type do not get a new bond in order to avoid creating multiple bonds bind_at_point_of_collision prevents sliding of the particles at the contact This is achieved by creating two virtual sites at the point of collision They are rigidly connectd to the colliding particles respectively A bond is then created between the virtual sites or an angular bond between the two real particles and the virtual particles In the latter case the virtual particles are the centers of the angle potentials particle 2 in the description of the angle potential see 5 5 Due to the rigid connection between each of the particles in the collision and its respective virtual site a sliding at the contact point is no longer possible See the documentation on rigid bodies for details In addition to the bond between the virtual sites the bond between the colliding particles is also created You can either use a real bonded interaction to prevent wobbling arou
101. anged either when the configure script is called with the option with myconfig see section 3 1 on page 20 or when make is called with the setting myconfig myconfig_header nET on pes 2 or mie The configuration header can be used to compile different binary versions of ESPResSo with a different set of features from the same source directory Suppose that you have a source directory srcdir and two build directories builddir1 and builddir2 that contain different configuration headers e builddir1 myconfig h define ELECTROSTATICS define LENNARD JONES e builddir2 myconfig h define LJCOS 24 Then you can simply compile two different versions of ESPResSo via cd builddir1 srcdir configure make cd builddir2 srcdir configure make 25 4 Setting up particles 4 1 part Creating single particles 4 1 1 Defining particle properties Syntax part pid pos x y z type typeid v vz vy vz fe fy fz bond bondid pid q charge quat q1 q2 q3 q4 omega x y z torque x y z rinertia x y 2 unjfix x y 2 ext_force x y 2 ext_torque z y z exclude pid2 exclude delete pid2 mass mass dipm moment dip dx dy dz virtual v vs_relative pid distance vs_auto_relate_to pid temp T gamma g rotation rot Required features ELECTROSTATICS ROTATION EXTERNAL_FORCES EXCLUSION 5mass DIPOLES VIRTUAL_SITES_cOM VIRTUAL_SITES_RELATIVE 9 LANGEVIN_PER_PARTICL
102. ariant 1 disables Coulomb interactions Variant 2 returns the current pa rameters of the coulomb interaction as a Tcl list using the same syntax as used to setup the method e g coulomb 1 0 p3m 7 75 8 5 0 1138 0 0 coulomb epsilon 0 1 n_interpol 32768 mesh_off 0 5 0 5 0 5 Variant 3 is the generic syntax to set up a specific method or its parameters the details of which are described in the following subsections Note that using the electrostatic interaction also requires assigning charges to the particles This is done using the part command to set the charge q e g inter coulomb 1 0 p3m tune accuracy le 4 part 0 q 1 0 part 1 q 1 0 5 7 1 Coulomb P3M Syntax inter coulomb lg p3m feu mesh mesh meshy meshz cao alpha Required features ELECTROSTATICS Description For this feature to work you need to have the fftw3 library installed on your system In ESPResSo you can check if it is compiled in by checking for the feature FFTW 61 This command activates the P3M method to compute the electrostatic interactions between charged particles The different parameters are described in more detail in 15 Teut The real space cutoff as a positive floating point number mesh The number of mesh points as a single positive integer mesh yz The number of mesh points in x y and z direction This is relevant for noncubic boxes cao The charge assignment order an integer between 0 and 7 alpha The Ewald parameter as
103. around 10 30 e ELECTROSTATICS This switches on the various electrostatics algorithms such as P3M See section 5 7 on page 61 for details on these algorithms e DIPOLES This activates the dipole moment property of particles In addition the various magnetostatics algorithms such as P3M are switched on See section 5 7 on page 61 for details on these algorithms e ROTATION Switch on rotational degrees of freedom for the particles as well as the corresponding quaternion integrator See section 4 1 1 on page 26 for details Note that when the feature is activated every particle has three additional degrees of freedom which for example means that the kinetic energy changes at constant temperature is twice as large e ROTATION_PER_PARTICLE Allows to set whether a particle has rotational degrees of freedom e LANGEVIN_PER_PARTICLE Allows to choose the Langevin temperature and friction coefficient per particle e ROTATIONAL_INERTIA e EXTERNAL_FORCES Allows to define an arbitrary constant force for each particle individually Also allows to fix individual coordinates of particles e g keep them at a fixed position or within a plane 197 e CONSTRAINTS Turns on various spatial constraints such as spherical compartments or walls This constraints interact with the particles through regular short ranged potentials such as the Lennard Jones potential See section 4 3 on page 36 for possible constraint forms
104. ary to determine the swap probablility are calculated like this proc swap id temp1 temp2 global config particle delete foreach p config id eval part p set epot expr analyze energy total analyze energy kinetic return expr epot temp1 expr epot temp2 J Note that only the potential energy is taken into account The temperature enters only indirectly through the inverse temperature prefactor see Eqn 7 1 The simulation is then started as follows One of the processes runs the command 93 for set T 0 T lt 3 set T expr T 0 01 lappend temperatures T parallel_tempering main load 30 values temperatures rounds 1000 init init swap swap perform perform This command turns the ESPResSo instance executing it into the master part of the parallel tempering simulation It waits until a sufficient number of clients has connected This are additional ESPResSo instances which are identical to the master script except that they execute parallel_tempering main connect host load 30 init init swap swap perform perform Here host is a variable containing the TCP IP hostname of the computer running the master process Note that the master process waits until enough processes have connected to start the simulation In the example there are 300 temperatures and each process including the master process will deal with 30 of them Therefore 1 master and 9 slave processes are required
105. ased analysis and core analysis is the tclcommand observable that allows use the return value of arbitrary Tcl functions also self written as input for the core analysis See more below 9 1 Observables 9 1 1 Introduction The first step of the core analysis is to tell ESPResSo to create an observable An observable in the sense of the core analysis can be considered as a rule how to compute a certain set of numbers from a given state of the system It does not refer to the numbers itself Creating an observable means just allocating the corresponding memory assigning a function to compute the observable value and reserving an id which will be used to refer to the observable In addition to the possibility to print the observable value return the observable value to the script interface the id of a core observable can be passed to another analysis function The observable value is computed from the current state of the system at the moment when it is needed e when requested explicitly by the user calling the observable print function or when requested automatically by some other analysis function Not all observables are implemented in parallel When performing a parallel compu tation too frequent updates to observables which are not implemented in parallel may produce a significant slowdown 115 9 1 2 Creating an observable To create a new observable use Syntax observable new name parameters Description U
106. aster In a non periodic system the DAWAANR method gives the exact result 5 8 4 Magnetic Dipolar Direct Sum MDDS Syntax inter magnetic lg mdds n_cut value_n_cut Required features DIPOLES MAGNETIC_DIPOLAR_DIRECT_SUM Description The command enables the magnetic dipolar direct sum The dipole dipole interaction is computed by explicitly summing over all pairs If the system is periodic in one or more directions the interactions with further value_n_cut replicas of the system in all periodic directions is explicitly computed 70 As it is very slow this method is not intended to do simulations but rather to check the results you get from more efficient methods like P3M 5 9 Special interaction commands 5 9 1 Tunable slip boundary interaction Syntax inter typel type2 tunable_slip T YL Teut t Ur Vy vz Required features TUNABLE_SLIP Description Simulating microchannel flow phenomena like the Plane Poiseuille and the Plane Couette Flow require accurate boundary conditions There are two main boundary conditions in use 1 slip boundary condition which means that the flow velocity at the the hydrody namic boundaries is zero 2 partial slip boundary condition which means that the flow velocity at the hydro dynamic boundaries does not vanish In recent years experiments have indicated that the no slip boundary condition is indeed usually not valid on the micrometer scale Instead it has to be replaced b
107. ate graphics primitives in VMD that represent some of the spatial constraints sphere rhomboid and cylinder at present When wait is provided ESPResSo will wait for wait seconds for VMD to connect and only continue after this time has passed or VMD has connected If start is 1 the default it will automatically try to start VMD and connect to the ESPResSo simulation otherwise it writes a corresponding script to the file filename vmd_start script that can be executed via VMD either from the command line vmd e vmd_start script or from the Tcl console of VMD with the command play vmd_start script 137 10 7 Errorhandling Errors in the parameters are detected as early as possible and hopefully self explanatory error messages returned without any changes to the data in the internal data of ESPResSo This include errors such as setting nonexistent properties of par ticles or simply misspelled commands These errors are returned as standard Tcl errors and can be caught on the Tcl level via catch script err When run noninteractively Tcl will return a nice stack backtrace which allows to quickly find the line causing the error However some errors can only be detected after changing the internal structures so that ESPResSo is left in a state such that integration is not possible without massive fixes by the users Especially errors occuring on nodes other than the primary node fall under this condition for example a broken
108. ated by direct summation of the near formula This means that in addition to the algorithm above one has to only a few things during the calculation of the particle and cell blocks and one additionally calculates the contributions of the image charges and puts them either in a separate array or for the boundary layers into two extra cell blocks outside the simulation box The entries in the separate array are then added up over all processors and stored in the terms of the lowest topmost layer This are all modifications necessary for the far formula part In addition to the far formula part there is an additional loop over the particles at the boundary to directly calculate their interactions with their images For details refer to Tyagi et al 55 E 3 MMM1D In one dimensionally periodic systems with z being the periodic coordinate the far formula looks like p 2 4uz Xpo Kolwp cos wz 2u log 2uzy F p z 8ru2 pro PA1 wp cos wz 2u F p 2 8rul Dogo PKo wp sin wz the near formula is 214 lp z uno E A a ot py 2uzY 2 Folp z u Enzo 2 Se O po F p z u es 3 A EA uzp yui z k z y Zkz k 1 C l r J where p denotes the xy distance of the particles As for the two dimensional periodic case the obtained energy is equal to the one dimensional Ewald sum Algorithmically MMMID is uninteresting since neither the near nor far formula allow
109. autocorrelation is to be calculated Calculate an autocorrelation function on a set of data mbtools utils distance posi pos2 e posi A position vector e pos2 A position vector Calculate the distance between two points whose position vectors are given mbtools utils distance_min posi pos2 e posi A position vector e pos2 A position vector Calculate the minimum image distance between two position vectors mbtools utils min_vec posi pos2 e pos A position vector e pos2 A position vector Calculate the minimum image vector from position vector2 to postition 1 i e posl pos2 mbtools utils normalize vec e vec The vector to be normalised Normalize a vector mbtools utils scalevec vec scale e vec The vector to be scaled e scale Scaling factor Multiply all elements of a vector by a scaling factor mbtools utils uniquelist original e original A list possibly containing duplicate elements Construct a list of all the unique elements in the original list removing all duplication 178 14 6 5 Miscellaneous procs mbtools utils trap_mols molstotrap e molstotrap A list of trap values for molecules This list would typically be obtained by calling mbtools get_trappedmols immediately after the system has been setup Set the trap value for a list of molecules mbtools utils isoutside pos zone e pos The point whose status is to be determined e zone This will be a tcl list The fi
110. bond or illegal parameter combinations For error conditions such as the examples given above a Tcl error message of the form tcl_error background 0 error_a error_b 1 error_c is returned Following possibly a normal Tcl error message after the background key word all severe errors are listed node by node preceeded by the node number A special error is lt consent gt which means that one of the slave nodes found exactly the same er rors as the master node This happens mainly during the initialization of the integrate e g if the time step is not set In this case the error message will be background_errors 0 time_step not set 1 lt consent gt In each case the current action was not fulfilled and possibly other parts of the internal data also had to be changed to allow ESPResSo to continue so you should really know what you do if you try and catch these errors 138 11 Auxilliary commands 11 1 Finding particles and bonds 11 1 1 countBonds Syntax countBonds particle ist Description Returns a Tcl list of the complete topology described by particle_list which must have the same format as the output of the command part see section 4 1 on page 26 The output list contains only the particle id and the corresponding bonding in formation thus it looks like e g 106 0 107 107 0 106 0 108 108 0 107 0 109 210 0 209 0 211 211 0 210 212 213 for a single chain of 106 monomers b
111. bose output time_vs_oop Calculates the orientational order parameter S for each lipid through a call to the espresso command analyze lipid_orient_order stress_tensor verbose output time_vs_stress_tensor Calculates all 9 elements of the pressure tensor for the system through a call to the espresso command analyze stress_tensor pressure verbose output time_vs_pressure Calculates the isotropic pressure through a call to analyze pressure Results are printed as a list of the various contributions in the following order p_inst total ideal FENE harmonic nonbonded Where p_inst is the instantaneous pressure obtained directly from the barostat stray verbose output time_vs_stray Calculates the number of stray lipids based on a call to analyze get_lipid_orients 14 4 3 Adding a new routine To add a new analysis routine you should create a new file called myanalysis tcl which will contain all of your code At the top of this file you should declare a namespace for your analysis code and include all of the internal variables inside that namespace as follows 165 namespace eval mbtools analysis myanalysis variable av_myresult variable av_myresult_i variable f_tvsresult variable verbose namespace export setup_myanalysis namespace export analyze_myanalysis namespace export printav_myanalysis namespace export resetav_myanalisis Import your new file into the analysis package by adding a
112. ce subtracted Lennard Jones tabulated tabulated bond Tunable slip boundary interactions volume conservation interactive mode internal distances within a chain invalidate_system Tcl command 90 label DP Dthermostat lb Tcl command 148 length unit 9 Lennard Jones cosine interaction Lennard Jones interaction local area conservation local stress tensor local_box_1 global variable Maggs method Magnetostatic interactions make max_cut global variable max_cut_bonded global variable max_cut_nonbonded global variable max_num_cells global variable Maxwell Equation Molecular Dynamics MDDS method MEMD metadynamics Tcl command min_num_cells global variable minimal particle distance MMMID method MMM2D method moment of inertia matrix momentum exchange method 226 Morse interaction MPI Multiple tau correlator 124 myconfig h n_part_types global variable NEMD nemd Tcl command node_grid global variable Non bonded interactions npt_p_ext global variable npt_p_inst global variable nptiso_gamma0 global variable nptiso_gammav global variable observable Tcl command Observables P3M method parallel_tempering Tcl command 91 part Tcl command 26 particle distance particle distribution particles in the neighbourhood pearl necklace structures periodicity global variable physical units 9
113. ce but this may strongly depend on the observables you are correlating For more information we recommend to read Ref or to perform your own tests compress1 and compress2 Are functions used to compress the data when going to the next level of the multiple tau correlator Different compression functions for different observables can be specified if desired otherwise the same function is used for both Default is discard which takes one of the observable values and discards the other one This is safe for all observables but produces poor statistics in the tail For some observables linear compression can be used which makes an average of two neighbouring values but produces systematic errors Depending on the observable the systematic error can be anything between harmless and disastrous For more information we recommend to read Ref or to perform your own tests 9 2 3 Inquiring about already existing correlations Syntax 1 correlation 2 correlation n_corr Description Variant 1 returns a tcl list of the defined correlations including their parameters Variant 2 returns the number of currently defined correlations 122 9 2 4 Collecting time series data for the correlation Syntax 1 correlation id autoupdate start stop 2 correlation id update 3 correlation id finalize Description Variant 1 is the recommended way of updating the correlations By specifying start or stop it starts or stops automatic
114. ce will help the algorithmic error as we can tune the speed of light to a higher value At the same time it increases the interpolation error at an even higher rate Therefore momentarily it is advisable to choose the lattice with a rather fine mesh of the size of the particles As a rule of thumb the error will then be less than 107 for the particle force For a more detailed description of the algorithm see appendix D on page 204 or the publications 33 38 5 7 6 Electrostatic Layer Correction ELC Please cite 6 BIBTEX key elc in file doc ug citations bib when using ELC and in addition BIBTEX key icelc in file doc ug citations bib if you use dielectric interfaces Syntax inter coulomb elc mazimal_pairwise_error gap_size far_cutoff noneutralization dielectric m dielectric contrasts A Aj Required features ELECTROSTATICS Description This is a special procedure that converts a 3d method to a 2d method in computa tional order N Currently it only supports P3M This means that you will first have to set up the P3M algorithm via inter coulomb p3m params before using ELC The algorithm is definitely faster than MMM2D for larger numbers of particles gt 400 at reasonable accuracy requirements The maximal pairwise error maximal_pairwise_error sets the LUB error of the force between any two charges without prefactors see the pa pers The algorithm tries to find parameters to meet t
115. chain e constraints If this option is specified the particle setup up tries to obey pre viously defined constraints see section 4 3 on page 36 4 2 2 counterions Setting up counterions Syntax counterions N start pid mode SAW RW shield trymax charge val type typeid Required features ELECTROSTATICS Description This command will create N counterions in the simulation box Arguments e start pid Sets the particle id of the first counterion It defaults to the cur rent number of particles i e counterions are placed after all previously defined particles e mode SAW RW shield trymax Specifies the setup method to place the counterions It defaults to SAW See the polymer command for a detailed descrip tion e charge val Specifies the charge of the counterions If not set it defaults to 1 0 e type typeid Specifies the particle type of the counterions It defaults to 2 4 2 3 salt Setting up salt ions Syntax salt Ny N_ start pid mode SAW RW shield trymax charges val val_ types typeid typeid_ rad r Required features ELECTROSTATICS Description Create N positively and N_ negatively charged salt ions of charge val and val_ within the simulation box Arguments e start pid Sets the particle id of the first positively charged salt ion It defaults to the current number of particles e mode SAW RW shield trymax Specifi
116. conditions Syntax lbboundary shape shape_args velocity vz vy vz Required features LB_BOUNDARIES Description If nothing else is specified periodic boundary conditions are assumed for the LB fluid The lbboundary command allows to set up other internal or external boundaries The lbboundary command syntax is very close to the constraint syntax as usually one wants the hydrodynamic boundary conditions to be shaped similarily to the MD boundaries Currently the shapes mentioned above are available and their syntax exactly follows the syntax of the constraint command For example lbboundary wall dist 1 5 normal 1 0 0 creates a planar boundary condition at distance 1 5 from the origin of the coordinate system where the half space x gt 1 5 is treated as normal LB fluid and the other half space is filled with boundary nodes Intersecting boundaries are in principle possible but must be treated with care In the current only partly satisfactory all nodes that are within at least one boundary are treated as boundary nodes Improving this is nontrivial and suggestions are very welcome Currently only the so called link bounce back algorithm for wall nodes is available This creates a boundary that is located approximately midway between the lattice nodes so in the above example this corresponds indeed to a boundary at x 1 5 Note that the location of the boundary is unfortunately not independent of the viscosity This can
117. ctions is presented in Section 157 14 External package mbtools mbtools is a set of tcl packages for setting up analyzing and running simulations of lipid membrane systems mbtools comes with a basic set of tabulated forces and potentials for lipid interactions and some example scripts to help explain the syntax of the commands If you make use of mbtools or of these potentials please acknowledge us with a citation to Cooke I R Kremer K and Deserno M 2005 Tunable generic model for fluid bilayer membranes Phys Rev E 72 011506 14 1 Introduction mbtools is located in the folder Espresso packages mbtools One of the main features of mbtools is the ability to easily create initial lipid config urations with interesting geometries These include flat membranes cylinders spheres toroids and randomly distributed gases Each of these shapes is referred to as a geom etry and any number of geometries can be combined in a single simulation Once the geometry has been chosen the user specifies the molecules which should be placed in this geometry For example one could choose sphere as a geometry and then define two different lipids and or a protein to be placed on the sphere Within reason e g size restrictions it should be possible to use any mixture of known molecule types on any geometry The molecule types available at present include proteins lipids of any length and spherical colloids mbtools includes
118. cules of each molecule type mbtools utils minpartid topo 176 e topo A valid topology Minimum particle id for the given topology mbtools utils minmoltype topo e topo A valid topology Minimum molecule type id for this topology mbtools utils listmoltypes topo e topo A valid topology Make a list of all the molecule types in a topology Makes a check for duplication which would occur for an unsorted topology mbtools utils listmollengths topo e topo A valid topology Works out the length number of atoms of each molecule type and returns a list of these lengths 14 6 4 Math procs mbtools utils dot_product AB Returns A dot B mbtools utils matrix_vec_multiply A B Return the product of a matrix A with a vector B mbtools utils calc_proportions ilist Calculate the number of times each integer occurs in the list ilist mbtools utils average data from to e data A list of numbers to be averaged e from Optional starting index in data e to Optional ending index in data Calculate the mean of a list of numbers starting from from going up to to mbtools utils stdev data from to 177 e data A list of numbers to find the std deviation of e from Optional starting index in data e to Optional ending index in data Calculate the standard deviation of a list of numbers starting from from going up to to mbtools utils acorr data e data Data for which an
119. curacy occurs at the compression step In principle one can use any value of m and p to tune the algorithm performance However it turns out that using a high m dilutes the data at high Therefore m 2 is hard coded in the ESPResSo correlator and cannot be modified by user The value of p remains an adjustable parameter which can be modified by user by setting tau_lin when defining a correlation In general one should choose p gt m to avoid loss of statistical accuracy Choosing p 16 seems to be safe but it may depend on the properties of the analyzed correlation functions A detailed analysis has been performed in Ref 42 The choice of the compression function also influences the statistical accuracy and can even lead to systematic errors The default compression function is discard2 which discards the second for the compressed values and pushes the first one to the higher level This is robust and can be applied universally to any combination of observables and correlation operation On the other hand it reduces the statistical accuracy as the compression level increases In many cases the average compression operation can be applied which averages the two neighbouring values and the average then enters the higher level preserving almost the full statistical accuracy of the original data In general if averaging can be safely used or not depends on the properties of the difference 125 1 1 1 zi Bi p Aisi O Bi p 1 zi
120. cutoff 71 radius of the constraint Lennard Jones interactions Otherwise there is no possibility that the particles feel the viscous layer Please read reference 47 before using this interaction 5 9 2 DPD interaction Syntax inter typel type2 inter_dpd gamma r_cut wf tgamma tr_cut twf Required features INTER_DPD Description This is a special interaction that is to be used in conjunction with the Dissipative Particle Dynamics algorithm 6 2 3J when the INTER_DPD implementation is used The parameters correspond to the parameters of the DPD thermostat but can be set individually for the different interactions 5 9 3 Fixing the center of mass Syntax inter typeid1 typeid1 comfixed flag Required features COMFIXED Description This interaction type applies a constraint on particles of type typeid1 such that during the integration the center of mass of these particles is fixed This is accomplished as follows The sum of all the forces acting on particles of type typeid1 are calculated These include all the forces due to other interaction types and also the thermostat Next a force equal in magnitude but in the oppositte direction is applied on the particles This force is divided equally on all the particles of type typeid1 since currently there is no mass concept in ESPResSo Note that the syntax of the declaration of comfixed interaction requires the same particle type to be input twice If different particle types are giv
121. d inter forcecap see section or on an individual level using the reap variable When reap is set and inter forcecap individual has been issued before the maximal force that is generated by this potential is the force at reap By default force capping is off i e the cap radius is set to 0 5 1 4 Lennard Jones cosine interaction Syntax 1 inter typel type lj cos O Teut Toff 2 inter typel type lj cos2 O Tog w Required features 2 Licos 2 Licos2 45 Description specifies a Lennard Jones interaction with cosine tail 49 between particles of the types typel and type2 The first variant behaves as follows Until the minimum of the Lennard Jones potential at Tmin rof 250 it behaves identical to the unshifted Lennard Jones potential chit 0 Between fmin and reut a cosine is used to smoothly connect the potential to 0 i e Vir Le cos a r rom 6 1 5 3 where a 7 Create Tot Tmin rot and 8 T Tmin Tot Q 1 In the second variant the cutoff radius is Teut Tmin w Where fmin rof 260 as in the first variant The potential between fmin and Toys is given by V r ecos e Tmin 5 4 For r lt fmin V r is implemented as normal Lennard Jones potential see equation 5 1 with Cshift 0 Only the second variant allows capping the force using inter forcecap see sec tion 9 9 5 5 1 5 Smooth step interaction Syntax inter typel type2 smoot
122. demo scripts if you want to see all frames of the simulation Variant 4 will terminate the IMD session This is normally not only nice but also the operating system will not free the port for some time so that without disconnecting for some 10 seconds you will not be able to reuse the port 10 6 2 Using IMD in VMD The PDB PSF files created by ESPResSo via the command writepsf and writepdb can be loaded into VMD This should bring up an initial configuration Then you can use the VMD console to execute the command imd connect host port where host is the host running the simulation and port is the port it listens to Note that VMD crashes if you do that without loading a structure file before For more information on how to use VMD to extract more information or hide parts of configuration see the VMD Quick Help 10 6 3 Automatically setting up a VMD connection Syntax prepare_vmd_connection filename wait start constraints Description To reduce the effort involved in setting up the IMD connection starting VMD and load ing the structure file ESPResSo provides the command prepare_vmd_connection It writes out the required vsf structure description file to filename vsf default for filename is vmd doing some nice stuff such as coloring the molecules bonds and counterions appropriately rotating your viewpoint and connecting your system to the visualization server If constraints is not 0 then the command will cre
123. dent of the bulk viscosity up to a scaling factor It is however known that the values of the viscosity does affect the quality of the implemented link bounce back method gamma_odd and gamma_even are the relaxation parameters for the kinetic modes Due to their somewhat obscure nature they are to be given directly in LB units Before running a simulation at least the following parameters must be set up agrid dens visc tau friction For the other parameters the following are taken bulk_viscosity 0 gamma_odd 0 gamma_even 0 fr fy fz 0 Syntax lbfluid print_interpolated_velocity xz y z Description This variant returns the velocity at point in countinous space This can make it easier to calculate flow profiles independent of the lattice constant Syntax lbfluid save_ascii_checkpoint filename 1bfluid save_binary_checkpoint filename lbfluid load_ascii_checkpoint filename lbfluid load_ binary_checkpoint filename Description The first two save commands save all of the LB fluid nodes populations to filename in ascii or binary format respectively The two load commands load the populations from filename This is useful for restarting a simulation either on the same machine or a different machine Some care should be taken when using the binary format as the format of doubles can depend on both the computer being used as well as the compiler 149 This is currently only implemented for the cpu version of LB 12
124. dial distribution function rdf of particles with types specified in part_type_list_a around particles with types specified in part_type_list_b The range is given by rmin and rmaz and is divided into rbins equidistant bins Output format The output corresponds to the blockfile format see section 10 2 on page 127 parameters r raf r y 8 1 12 Structure factor Syntax analyze structurefactor type order Description Returns the spherically averaged structure factor S q for particles of a given type type The S q is calculated for all possible wave vectors Tr lt q lt order Do not chose parameter order too large because the number of calculations grows as order Output format The output corresponds to the blockfile format see section 10 2 on page 127 q value S q _value y 100 8 1 13 Van Hove autocorrelation function G r t Syntax analyze vanhove type rmin rmax rbins tmaz Description Returns the van Hove auto correlation function G r t and the mean square displacement msd t for particles of type ptype for the configurations stored in the array configs This tool assumes that the configurations stored with analyze append see section are stored at equidistant time intervals G r t is calculated for each multiple of this time intervals For each time t the distribution of particle displacements is calculated according to the specification given by rmin rmax and rbins Optional arg
125. ds some more elaborate organisation which will be presented here A particle itself is represented by a structure Particle consisting of several substructures e g ParticlePosition ParticleForce or ParticleProperties which in turn represent basic physical properties such as position force or charge The particles are organised in one or more particle lists on each node called Cell cells The cells are arranged by several possible systems the cellsystems as described above A cell system defines a way the particles are stored in ESPResSo i e how they are distributed onto the processor nodes and how they are organised on each of them Moreover a cell system also defines procedures to efficiently calculate the force energy and pressure for the short ranged interactions since these can be heavily optimised depending on the cell system For example the domain decomposition cellsystem allows an order N interactions evaluation Technically a cell is organised as a dynamically growing array not as a list This ensures that the data of all particles in a cell is stored contiguously in the memory The particle data is accessed transparently through a set of methods common to all cell systems which allocate the cells add new particles retrieve particle information and are responsible for communicating the particle data between the nodes Therefore most portions of the code can access the particle data safely without direct knowledge of the current
126. duced_charges epsilons icc_epsilons normals icc_normals areas icc_areas sigmas icc_sigmas 5 8 Dipolar interaction Syntax 1 inter magnetic 0 0 2 inter magnetic 3 inter magnetic parameters Description These commands can be used to set up magnetostatic interactions which is defined as follows dd B PEP UP P3M 5 loko 7 j i j 5 37 where r r lg is a dimensionless parameter similar to the Bjerrum length in electrostatics which helps to tune the effect of the medium on the magnetic interaction between two magnetic dipoles Computing magnetostatic interactions is computationally very expensive ESPResSo features some state of the art algorithms to deal with these interactions as efficiently as possible but almost all of them require some knowledge to use them properly Unedu cated use can result in completely unphysical simulations The commands above work as their couterparts for the electrostatic interactions see section 5 7 1 on page 61 Variant 1 disables dipolar interactions Variant 68 2 returns the current parameters of the dipolar interaction as a Tcl list using the same syntax as used to setup the method e g coulomb 1 0 p3m 7 75 8 5 0 1138 0 0 coulomb epsilon 0 1 n_interpol 32768 mesh_off 0 5 0 5 0 5 Variant 3 is the generic syntax to set up a specific method or its parameters the details of which are described in the following subsections Note that using the magnetosta
127. dynamics algorithm J Phys Condens Matter 16 38 3999 4020 September 2004 D 220 39 40 41 42 43 44 45 46 47 48 49 50 51 S Poblete M Praprotnik K Kremer and Luigi Delle Site Coupling different levels of resolution in molecular simulations Jour Chem Phys 132 11 114101 2010 Matej Praprotnik Luigi Delle Site and Kurt Kremer Adaptive resolution molecular dynamics simulation Changing the degrees of freedom on the fly The Journal of Chemical Physics 123 22 224106 14 2005 6 5 Matej Praprotnik Luigi Delle Site and Kurt Kremer Multiscale simulation of soft matter From scale bridging to adaptive resolution Annual Review of Physical Chemistry 59 1 545 571 2008 6 5 Jorge Ramirez Sathish K Sukumaran Bart Vorselaars and Alexei E Likhtman Efficient on the fly calculation of time correlation functions in computer simulations J Chem Phys 133 15 154103 OCT 21 2010 ISSN 0021 9606 doi 10 1063 1 3491098 D Roehm and A Arnold Lattice boltzmann simulations on GPUs with ESPResSo Eur Phys J ST 210 73 88 2012 Michael Rubinstein and Ralph H Colby Polymer Physics Oxford University Press Oxford UK 2003 K Schatzel M Drewel and S Stimac Photon correlation measurements at large lag times improving statistical accuracy Journal of Modern Optics 35 4 711 718 APR 1988 ISSN 0950 0340 doi 10 1080 09500348814550731
128. e TUNABLE_SLIP Switch on tunable slip conditions for planar wall boundary condi tions See section 5 9 1 on page 71 for details e MASS Allows particles to have individual masses Note that some analysis proce dures have not yet been adapted to take the masses into account correctly e EXCLUSIONS Allows to exclude specific short ranged interactions within molecules e COMFORCE Allows to pull apart groups of particles e COMFIXED Allows to fix the center of mass of all particles of a certain type e MOLFORCES e BOND_CONSTRAINT Turns on the RATTLE integrator which allows for fixed lengths bonds between particles e VIRTUAL_SITES_COM Virtual sites are particles the position and velocity of which is not obtained by integrating equations of motion Rather they are placed using the position and orientation of other particles The feature VIRTUAL_SITES_COM allows to place a virtual particle into the center of mass of a set of other particles See section for details e VIRTUAL_SITES_RELATIVE Virtual sites are particles the position and velocity of which is not obtained by integrating equations of motion Rather they are placed using the position and orientation of other particles The feature VIRTUAL_SITES_RELATIVE allows for rigid arrangements of particles See section 4 4 for details e VIRTUAL_SITES_NO_VELOCITY e VIRTUAL_SITES_THERMOSTAT e THERMOSTAT_IGNORE_NON_VIRTUAL e BOND_VIRTUAL e MODES e ADRESS e METADYNAMICS e LA
129. e center of a molecule perform the following steps in that order 1 Create a particle of the desired type for each molecule It should be placed at least roughly in the center of the molecule to make sure it s on the same node as the other particles forming the molecule in a simulation with more than one cpu 2 Make it a virtual site using art pid virtual 1 P p 3 Declare the list of molecules and the particles they consist of analyze set molid listofparticleids The lists of particles in a molecule comprise the non virtual particles and the vir tual site 4 Assign to all particles that belong to the same molecule a common molecule id art pid mol molid P 5 Update the position of all virtual particles optional integrate 0 39 4 4 2 Rigid arrangements of particles The relative implementation of virtual sites allows for the simulation of rigid arrange ments of particles It can be used e g for extended dipoles and raspberry particles but also for more complex configurations Position and velocity of a virtual site are obtained from the position and orientation of exactly one non virtual particle which has to be placed in the center of mass of the rigid body Several virtual sites can be related to one and the same non virtual particle The position of the virtual site is given by E E On O E Jd 4 1 where zx is the position of the non virtual particle On is the orientation of the
130. e contained within path to be preceded by a prefix and having postfix before the suffix e g timeStamp scripts config gz DH863 001 gz returns scripts DH863_config001 gz If postfix is 1 the current date is used in the format y m 4d This would results in scripts DH863_config021022 gz on October 22nd 2002 11 2 Additional Tcl math functions The following procedures are found in scripts ABHmath tcl e CONSTANTS PI returns 7 with 16 digits precision KBOLTZ Returns Boltzmann constant in Joule Kelvin ECHARGE Returns elementary charge in Coulomb NAVOGADRO Returns Avogadro number SPEEDOFLIGHT Returns speed of light in meter second EPSILONO 140 Returns dielectric constant of vaccum in Coulomb2 Joule meter ATOMICMASS Returns the atomic mass unit u in kilogramms e MATHEMATICAL FUNCTIONS sqr lt arg gt returns the square of arg min lt argi gt lt arg2 gt returns the minimum of arg1 and arg2 max lt argl gt lt arg2 gt returns the maximum of arg1 and arg2 sign lt arg gt returns the signum function of arg namely 1 for arg gt 0 1 for lt 0 and 0 otherwise e RANDOM FUNCTIONS gauss_random returns random numbers which have a Gaussian distribution dist_random lt dist gt max returns random numbers in the interval 0 1 which have a distribution ac cording to the distribution function p x dist which has to be given as a tcl list containing equally spaced
131. e default behaviour is to print from the main mbtools namespaces and the global namespace mbtools system_generation mbtools utils mbtools analysis Note that children of these namespaces must be explicitly enabled All message types except debug are also enabled by default The following commands allow this default behaviour to be changed immsg setnamespaces namespacelist e namespacelist A list of all namespaces from which messages are to be printed Allows control over which namespaces messages can be printed from mmsg enable type e type A string indicating a single message type to enable Allowable values are err debug send and warn Allows particular message types to be enabled For example one could enable debug output with mmsg enable debug mmsg disable type e type A string indicating a single message type to disable Allowable values are err debug send and warn Allows particular message types to be disabled For example one could disable warning output with mmsg enable warn 181 15 Under the hood e Implementation issues that are interesting for the user e Main loop in pseudo code for comparison 15 1 Internal particle organization Since basically all major parts of the main MD integration have to access the particle data efficient access to the particle data is crucial for a fast MD code Therefore the particle data nee
132. e equivalent part O bond 44 1 2 part O bond 44 2 1 part 1 bond 44 0 2 5 4 4 Global area conservation Syntax inter bondid area_force_global 9 kag Description This type of interaction is available solely for closed 3D immersed objects 57 The conservation of local area is sometimes too restrictive Denote by S the current surface of the immersed object by Sg the surface in the relaxed state and define AS S Sy The global area conservation force is defined as AS Here the above mentioned force divided by 3 is added to all three particles A Cc Again this interaction is symmetric as is the area_force_local 5 4 5 Volume conservation Syntax inter bondid volume_force V ky Description This type of interaction is available solely for closed 3D immersed objects The deviation of the objects volume V is computed from the volume in the resting shape AV V V For each triangle the following force is computed F ABC ko hr Sano NABC 5 26 where SABc is the area of triangle ABC nac is the normal unit vector of plane ABC and ky is the volume constraint coefficient The volume of one immersed object is computed from V Y Sago nasc hasc 5 27 ABC where the sum is computed over all triangles of the mesh and hago is the normal vector from the centroid of triangle ABC to any plane which does not cross the cell The force F ABC is equally distributed to all three vertices A B C This i
133. e explicit zone For more details on the method itself see 40 41 39 And for further information about the technical implementation see 27 6 6 CUDA Syntax 1 cuda list 2 cuda setdevice id cuda getdevice 3 g Description This command can be used to choose the GPU for all subsequent GPU computations Note that due to driver limitations the GPU cannot be changed anymore after the first GPU using command has been issued for example 1bfluid If you do not choose the GPU manually before that CUDA internally chooses one which is normally the most powerful GPU available but load independent Variant 1 lists the available devices by their ids and brand names Variant 2 allows to choose the device by its id which can be determined using cuda list or for example the deviceQuery example code in the CUDA SDK Variant 3 finally gives the id of the currently active GPU 6 7 Creating bonds when particles collide Please cite 7 BIBTEX key espresso2 in file doc ug citations bib when using dynamic bonding 83 on_collision on_collision off 2 4 on_collision exception bind_centers d bond1 5 on_collision exception bind_at_point_of_collision d bondi bond2 type Description With the help of the feature COLLISION_DETECTION bonds between particles can be created automatically during the simulation every time two particles collide This is useful for simulations of chemical reactions
134. e inverted speed of light e mesh is the number of mesh points for the interpolation of the electromagnetic field in one dimension eso is the background dielectric permittivity at infinity This defaults to metallic boundary conditions to match the results of P3M The arising self interactions are treated with a modified version of the exact solution of the lattice Green s function for the problem Currently forces have large errors for two particles within the same lattice cube This may be fixed in future development but right now leads to the following rule of thumb for the parameter choices e The lattice should be of the size of your particle size i e the lennard jones epsilon That means mesh box_l lj_sigma e The integration timestep should be in a range where no particle moves more than one lattice box i e lennard jones sigma per timestep e The speed of light should satisfy the stability criterion c lt a dt where a is the lattice spacing and dt is the timestep For the second parameter this means f mass gt dt a 65 The main error of the MEMD algorithm stems from the lattice interpolation and is proportional to the lattice size in three dimensions which means Ajattice X Q Without derivation here the algorithmis error is proportional to 1 c where c is the adjustable speed of light From the stability criterion this yields A maggs As a B dt a 5 36 This means that increasing the latti
135. e molecule at the desired position with the desired orientation 171 14 5 3 Adding a new geometry To create a routine for setting up a system with a new type of geometry mygeom Start by creating a new file mygeom tcl inside the system_generation directory The new file should declare a new namespace mygeom as a sub namespace of mbtools system_ generation and export the proceedure create_mygeom Thus your mygeom tcl file should begin with the lines namespace eval mbtools system_generation mygeom namespace export create_mygeom Import your new file into the system_generation package by adding a line like the following to the system_generation tcl file source file join file dirname info script mygeom tcl You then need to implement the create_mygeom proceedure within your new names pace as follows mbtools system_generation mygeom create_mygeom args 14 5 4 Available molecule types lipid typeinfo moltypeid lipid particletypelist bondtypelist e particletypelist A list of the particle types for each atom in the lipid The particles are placed in the order in which they appear in this list e bondtypelist A list of two bondtypeids The first id is used for bonds between consecutive beads in the lipid The second bondtypeid defines the pseudo bending potential which is a two body bond acting across beads separated by exactly one bead Places atoms in a line to create a lipid molecule ho
136. e new values respectively re initialising the random num ber generators on each node Note that this is automatically done on invoking Espresso however due to that your simulation will always start with the same ran dom sequence on any node unless you use this tcl command to reset the sequences seeds 144 e Since internally the random number generators random sequences are not based on mere seeds but rather on whole random number tables to recover the exact state of the random number generators at a given time during the simulation run e g for saving a checkpoint requires knowledge of all these values They can be accessed by t_random stat which returns a tcl list with all status informations for any node e g 8 nodes gt approx 350 parameters To overwrite those internally in Espresso e g upon restoring a checkpoint submit the whole list back using t_random stat lt status list gt with status list being the tcl list mentioned above without any braces Be careful A complete recovery of the current state of the simulation is only possible if you make sure to include a call to The invalidate_system command after you saved the checkpoint tcl_checkpoint_set will do this automatically for you because the integration algorithm re uses the old forces calculated in the previous time step if something has changed in the system or if it has just been read from a file the forces are re derived including applicati
137. e nsquared only with an odd number of nodes if with multiple processors at all 6 4 3 Layered cell system Syntax cellsystem layered n_layers Description This selects the layered cell system which is specifically designed for the needs of the MMM2D algorithm Basically it consists of a nsquared algorithm in x and y but a domain decomposition along z i e the system is cut into equally sized layers along the z axis The current implementation allows for the cpus to align only along the z axis therefore the processor grid has to have the form 1x1xN However each processor may be responsible for several layers which is determined by n_layers i e the system is split into N n_layers layers along the z axis Since in x and y direction there are no processor boundaries the implementation is basically just a stripped down version of the domain decomposition cellsystem 6 5 AdResS Please cite BIBTEX key adress in file doc ug citations bib when using AdRess 82 Syntax adress set topo kind width width hybrid_width center x R_x wf wf Required features ADRESSO Description where kind determines the type of AdResS simulation 0 disabled 1 constant weight function 2 one dimensional splitting 3 spherical splitting wf the type of weighting function 0 standard 1 user defined width and hybrid_width are the widths of the explicit and hybrid regions respectively and R_z is the x position of the center of th
138. e particle lists for each cell there is only one particle list per node and a the cells actually only contain pointers into this particle list This has the advantage that when particles are moved from one cell to another on the same processor only the pointers have to be updated which is much less data 4 rsp 8 bytes than the full particle structure around 192 bytes depending on the features compiled in The data storage scheme of ESPResSo however requires to always move the full particle data Nevertheless from our experience the second approach is 2 3 times faster than the classical one To understand this one has to know a little bit about the architecture of modern computers Most modern processors have a clock frequency above 1GHz and are able to execute nearly one instruction per clock tick In contrast to this the memory runs at a clock speed around 200MHz Modern double data rate DDR RAM transfers up to 3 2GB s at this clock speed at each edge of the clock signal 8 bytes are transferred But in addition to the data transfer speed DDR RAM has some latency for fetching the data which can be up to 50ns in the worst case Memory is organised internally in pages or rows of typically 8KB size The full 2 x 200 MHz data rate can only be achieved if the access is within the same memory page page hit otherwise some latency has to be added page miss The actual latency depends on some other aspects of the memory organisation which will
139. e very primitive nsquared cellsystem which calculates the interactions for all particle pairs Therefore it loops over all particles giving an unfavorable computation time scaling of N However algorithms like MMMID or the plain Coulomb interaction in the cell model require the calculation of all pair interactions In a multiple processor environment the nsquared cellsystem uses a simple particle balancing scheme to have a nearly equal number of particles per CPU i e n nodes have m particles and p n nodes have m 1 particles such that n m p n m 1 N the total number of particles Therefore the computational load should be balanced fairly equal among the nodes with one exception This code always uses one CPU for the interaction between two different nodes For an odd number of nodes this is fine because the total number of interactions to calculate is a multiple of the number of nodes but for an even number of nodes for each of the p 1 communication rounds one processor is idle E g for 2 processors there are 3 interactions 0 0 1 1 0 1 Naturally 0 0 and 1 1 are treated by processor O and 1 respectively But the 0 1 interaction is treated by node 1 alone so the workload for this node is twice as high For 3 processors the interactions are 0 0 1 1 2 2 0 1 1 2 0 2 Of these interactions node 0 treats 0 0 and 0 2 node 1 treats 1 1 and 0 1 and node 2 treats 2 2 and 1 2 Therefore it is highly recommended that you us
140. ed behaviour if you are sure not to violate the dissipation fluctuation theorem e DPD does not take into account any internal rotational degrees of freedom of the particles if ROTATION is switched on Up to the current version DPD only acts on the translatorial degrees of freedom 78 Transverse DPD thermostat This is an extension of the above standard DPD ther mostat 28 which dampens the degrees of freedom perpendicular on the axis between two particles To switch it on the feature TRANS_DPD is required instead of the feature DPD The dissipative force is calculated by FR Cw rig I Fag O Fag Dig The random force by Ej ow rag I fig O fig Oy The parameters tgamma tr_cut define the strength of the friction and the cutoff in the same way as above Note This thermostat does not conserve angular momentum Interaction DPD Syntax thermostat inter_dpd temperature Required features INTER_DPD Description Another way to use DPD is by using the interaction DPD In this case DPD is imple mented via a thermostat and corresponding interactions The above command will set the global temperature of the system while the friction and other parameters have to be set via the command inter inter_dpd see 5 9 2 on page 72 This allows to set the friction on a per interaction basis Other DPD extensions The features DPD_MASS_RED or DPD_MASS_LIN make the friction constant mass depen dent gt M
141. efer to the publications 38 The method is intimately related to the Car Parrinello approach while being equivalent to solving Maxwell s equations with freely adjustable speed of light D 1 Equations of motion Denoting the particle masses with m their charges with qi their coordinates and mo mentum with 7 and p respectively the interparticle potential of non electromagnetic type with U for the coupled system of charges and fields we write the following equa tions of motion 1 gt 23 D 1 T me D 1 2 OU E Pi TOF Gk Ti Pi fi D 2 Ti i A E D 3 Se 3 l gt E eV x Wx A D 4 0 where ey is the vacuum dielectric constant c the speed of light A the vector potential E the electric field j the current density is the particle friction constant and f is a random force satisfying the standard fluctuation dissipation theorem SOL 0 kB ij5a95 t t D 5 where a and P denote Cartesian indices If we introduce the vector B V x A the system of equations can be rewritten in a form similar to the usual Maxwell equations Currently in ESPResSo the version with B and E is implemented 204 D 2 Discretization For implementation on the computer the equations need to be discretized with respect to both space and time We consider a domain of physical space as being an affine space and divide it into subdomains of contiguous cells of cubic shape The charges live on
142. elation time and Q denotes a quality mea sure i e the probability to find a x fit of the replica estimates The function returns an error message if the windowing failed or if the error in one of the replica is to large 114 9 Analysis in the core Analysis in the core is a new concept introduced in ESPResSo since version 3 1 It was motivated by the fact that sometimes it is desirable that the analysis functions do more than just return a value to the scripting interface For some observables it is desirable to be sampled every few integrations steps In addition it should be possible to pass the observable values to other functions which compute history dependent quantities such as correlation functions All this should be done without the need to interrupt the integration by passing the control to the script level and back which produces a significant overhead when performed too often Some observables in the core have their corresponding counterparts in the Tcl ob servables of the analyze command described in Chapter However only the core observables can be used on the fly with the toolbox of the correlator and on the fly analysis of time series Similarly some special cases of using the correlator have their redundant counterparts in the analysis in Tcl Chapter 8 but the correlator provides a general and versatile toolbox which can be used with any implemented core observables The only trick to bridge the gap between Tcl b
143. elling and simulation of processes in microfluidic devices for biomedical applications Computers amp Mathematics with Applications 64 3 278 288 2012 M Dao C T Lim and S Suresh Mechanics of the human red blood cell deformed by optical tweezers J Mech Phys Solids 51 2259 2280 2003 M Deserno Counterion condensation for rigid linear polyelectrolytes PhD thesis Universitat Mainz 2000 5 7 1 M Deserno and C Holm How to mesh up Ewald sums i J Chem Phys 109 7678 1998 5 71 E31 M Deserno and C Holm How to mesh up Ewald sums ii J Chem Phys 109 7694 1998 M Deserno C Holm and H J Limbach Molecular Dynamics on Parallel Com puters chapter How to mesh up Ewald sums World Scientific Singapore 2000 M Doi and S F Edwards The theory of polymer dynamics Oxford Science Publi cations 1986 M M Dupin I Halliday C M Care and L Alboul Modeling the flow of dense suspensions of deformable particles in three dimensions Phys Rev E Stat Nonlin Soft Matter Phys 75 066707 2007 P P Ewald Die berechnung optischer und elektrostatischer gitterpotentiale Ann Phys 64 253 287 1921 5 7 1 Daan Frenkel and Berend Smit Understanding Molecular Simulation Academic Press San Diego second edition 2002 Gary S Grest and Kurt Kremer Molecular dynamics simulation for polymers in the presence of a heat bath Phys Rev A 33 5 3628 31 1986 Owen A Hickey Christian Ho
144. en in the input the program exits with an error message flag can be set to 1 which turns on the interaction or 0 to turn off the interaction 5 9 4 Pulling particles apart Syntax inter typeid1 typeid2 comforce flag dir force fratio Required features COMFORCE Description The comforce interaction type enables one to pull away particle groups of two different types It is mainly designed for pulling experiments on bundles Within a bundle of molecules of type number typeid1 lets mark one molecule as of type typeid2 Using 72 comforce one can apply a force such that t2 can be pulled away from the bundle The comforcerlag is set to 1 to turn on the interaction and to 0 otherwise The pulling can be done in two different directions Either parallel to the major axis of the bundle dir 0 or perpendicular to the major axis of the bundle dir 1 force is used to set the magnitude of the force fratio is used to set the ratio of the force applied on particles of typeid1 vs typeid2 This is useful if one has to keep the total applied force on the bundle and on the target molecule the same A force of magnitude force is applied on typeid2 particles and a force of magnitude force fratio is applied on typeid1 particles 5 9 5 Capping the force during warmup Syntax 1 inter forcecap Fmax individual Description Non bonded interactions are often used to model the hard core repulsion between par ticles Most of
145. eo Hookean behaviour 0 5 2 5 _ AAB t AAB FA MBE 5 22 A The stretching force acts between two particles and is symmetric Therefore if an inter action is defined by inter 1 stretching_force 2 0 4 0 then the following two commands 50 part 42 bond 1 43 part 43 bond 1 42 are equivalent 5 4 2 Bending force Syntax inter bondid bending_force 0 ky Description The tendency of an elastic object to maintain the resting shape is governed by prescrib ing the prefered angles between the neighbouring triangles of the mesh This type of interaction is available for closed 3D immersed objects as well as for 2D sheet flowing in the 3D flow Denote by 6 the angle between two triangles in the resting shape For closed immersed objects you always have to set the inner angle The deviation of this angle A8 0 0 is computed and defines two bending forces for two triangles A BC and Ag BC A0 Fyi A BC ky Hy MAiBC 5 23 Here na Bc is the unit normal vector to the triangle A BC The force Fy A BC is assigned to the vertex not belonging to the common edge The opposite force divided by two is assigned to the two vertices lying on the common edge This procedure is done twice for i 1 and for i 2 Unlike the stretching force the bending force is strictly asymmetric After creating an interaction inter 33 bending_force 0 7 4 0 it is important how the bond is created Particles need to be
146. epare for vmd output 175 14 6 2 Warmup commands mbtools utils warmup steps times mindist arg cfgs arg outputdir arg vmdflag arg startcap arg capgoal arg e steps number of integration steps used in each call to integrate e times number of times to call the integrate function during warmup e mindist 0 Terminate the warmup when the minimum particle distance is greater than this criterion A value of 0 default results in this condition being ignored If a condition is imposed this routine can become very very slow for large systems e cfgs 1 Write out a configuration file every cfgs calls to integrate e outputdir The directory for writing output e umdflag offline If this flag is set to offline default pdb files will be generated for each configuration file generated e startcap 5 Starting value for the forcecap e capgoal 1000 For the purposes of calculating a cap increment this value is used as a goal The final forcecap will have this value Perform a series of integration steps while increasing forcecaps from an initially small value 14 6 3 Topology procs mbtools utils maxpartid topo e topo A valid topology Find the maximum particle id in a given topology mbtools utils maxmoltypeid topo e topo A valid topology Find the maximum molecule type id mbtools utils listnmols topo e topo A valid topology Construct a list with the number of mole
147. er id which has been assigned to it Its further arguments are described below Arguments e obsi and obs2 are ids of the observables A and B that are to correlated The ids have to refer to existing observables which have been previously defined by the observable command Some observables are already implemented and others can be easily added This can be done with very limited ESPResSo knowledge just by following the implementations that are already in If obs2 is omitted autocorrelation of obs1 is calculated by default e corr_operation The operation that is performed on A t and B t T to obtain C 7 The following operations are currently is available e scalar_product Scalar product of A and B i e C Y A B e componentwise_product Comnponentwise product of A and B i e Ci A B e square_distance_componentwise Each component of the correlation vector is the square of the difference between the corresponding components of the observables i e Ci A B Example when A is particle_positions it produces the mean square displacement for each component separately e complex_conjugate_product e fcs_acf WrWy Wz Fluorescence Correlation Spectroscopy FCS autocorrelation function i e car oo E gt aay where Az T ai 0 xi t is the square discplacement of particle i in the x direction and wz is the beam waist of the intensity profile of the exciting laser beam 9 2 9 3 2
148. er installations sometimes require np instead of n or mpirun instead of mpiexec 3 4 myconfig h Activating and deactivating features ESPResSo has a large number of features that can be compiled into the binary However it is not recommended to actually compile in all possible features as this will slow down ESPResSo significantly Instead compile in only the features that are actually required A strong gain in speed can be achieved by disabling all non bonded interactions except for a single one e g LENNARD_JONES For the developers it is also possible to turn on or off a number of debugging messages The features and debug messages can be controlled via a configuration header file that contains C preprocessor declarations Appendix B lists and describes all available features When no configuration header is provided by the user a default header found in src myconfig default h will be used that turns on the default features The file myconfig sample h in the source directory contains a list of all possible features that can be copied into your own configuration file When you distinguish between the build and the source directory the configuration header can be put in either of these Note however that when a configuration header is found in both directories the one in the build directory will be used By default the configuration header is called myconfig h The name of the configu ration header can be ch
149. ere in the order they appear in the topology 169 Creates a spherical cap which is part of a vesicle of a radius r by placing molecules in an ordered manner at uniform density on the surface of the sphere Molecules are assumed to have a uniform cross sectional area and closely matched though not identical lengths If the option half is defined the radius of the vesicle will depend on the number of lipids and the area per lipid torus C arg initarea arg ratio arg bondl arg shuffle e c 0 0 0 0 0 0 The location of the center of the torus relative to the center of the box e initarea 1 29 An initial guess for the area per lipid This guess is used to compute initial radii based on the number of lipids This initial guess is then iteratively refined until all lipids can be fit uniformly on the torus e ratio 1 4142 Ratio of major toroidal radius to minor toroidal radius Default value is for the Clifford torus e shuffle shuffle the topology prior to placing the lipids This is required for a random lipid distribution because otherwise the lipids will be placed on the torus in the order they appear in the topology Creates a toroidal vesicle by placing molecules in an ordered manner at uniform density on the surface of the torus Molecules are assumed to have a uniform cross sectional area and closely matched though not identical lengths The two radii of the torus will depend on the number of lipids the area per lipid
150. ergy obtained by the spherical summation limit The deeper reason for this is that in some sense the electrostatic sum is absolutely convergent 5 The near formula is used for particles with a small distance along the z axis for all other particles the far formula is used Below is shown that the far formula can be evaluated much more efficiently however its convergence breaks down for small z distance To efficiently implement MMM2D the layered cell system is required which splits up the system in equally sized gaps along the z axis The interaction of all particles in a layer S with all particles in the layers S 1 5 S 1 is calculated using the near formula for the particles in layers 1 S 2 and in layers S 2 N the far formula is used The implementation of the near formula is relatively straight forward and can be treated as any short ranged force is treated using the link cell algorithm here in the layered variant The special functions in the formula are somewhat demanding but for the polygamma functions Taylor series can be achieved which are implemented in mmm common h The Bessel functions are calculated using a Chebychev series The treatment of the far formula is algorithmically more complicated For a particle i in layer S the formula can product decomposed as in Djersses 1 005 EEE C08 wp i 23 cos ugly 45 y cos w 2 cos wqyi gt jels S lt S 1 qje Tra i Cos Wpx cos wyy Ga cos wpT sin wgy
151. es the setup method to place the counterions It defaults to SAW See the polymer command for a detailed descrip tion e charge val val_ Sets the charge of the positive salt ions to val and the one of the negatively charged salt ions to val_ If not set the values default to 1 0 and 1 0 respectively e type typeid typeid_ Specifies the particle type of the salt ions It defaults to 3 respectively 4 32 e rad r The salt ions are only placed in a sphere with radius r around the origin 4 2 4 diamond Setting up diamond polymer networks Syntax diamond a bond length monomers_per_chain counterions Nc charges Valnode Valmonomer vale distance denarged 1 nonet Required features 1 ELECTROSTATICS Description Creates a diamond shaped polymer network with 8 tetra functional nodes connected by 2 x8 polymer chains of length monomers_per_chain in a unit cell of length a Chain monomers are placed at a mutual distance bond_length along the vector connecting network nodes The polymer is created starting from particle ID 0 Nodes are assigned type 0 monomers both charged and uncharged are type 1 and counterions type 2 For inter particle bonds interaction 0 is taken which must be a two particle bond Figure 4 1 Diamond like polymer network with monomers_per_chain 15 Arguments ea Determines the size of the of the unit cell e bond_length Specifies the bond length of the polymer chains connecting the
152. esSo does not predefine any units While most MD programs specify a set of units like for example that all lengths are measured in ngstr m or nanometers times are measured in nano or picoseconds and energies are measured in kJ mol ESPResSo does not do so Instead the length time and energy scales can be freely chosen by the user A length of 1 0 can mean a nanometer an ngstr m or a kilometer depending on the physical system that the user has in mind when he writes his ESPResSo script The user can choose the unit system that suits the system best When creating particles that are intended to represent a specific type of atoms one will probably use a length scale of Angstrom This would mean that e g the parameter o of the Lennard Jones interaction between two atoms would be set to twice the van der Waals radius of the atom in ngstr m Alternatively one could set to 2 0 and measure all lengths in multiples of the van der Waals radius The second choice to be made is the energy and time scale One can for example choose to set the Lennard Jones parameter e to the energy in EL Then all energies will be measured in that unit Alternatively one can choose to set it to 1 0 and measure everything in multiples of the van der Waals binding energy As long as one remains within the same unit system throughout the whole ESPResSo script there should be no problems 1 5 Requirements The following libraries and tools are
153. estricts observable calculation to a given list of particle id s The id list is a tcl list of existing particle ids 9 1 7 Profile specifications Profiles are specified by giving the spacial area that is to be profiled and the number of bins in each spacial direction The area to be analyzed is characterized by minzx mazz miny mazry and minz maxz The defaults correspond to the box size when the observ able is created The bin size in each direction defaults to 1 and can be change with the parameter xbins ybins zbins Changing one two or three of them to a value gt 1 with thus create a one two or three dimensional map of the desired quantity The full syntax thus reads as Syntax observable new needs_profile_specs other parameters minx minz maxx mazz miny miny maxy mazy minz minz maxz mazz xbins zbins ybins ybins zbins zbins Description Radial profiles allow to do the same as usual profiles except the coordinate system is a cylindrical one and the binning is done in the cylindrical coordinates defined with the axis in z direction This is very helpful if the symmetry of the system is cylindrical The spacial are is characterized by a center default to the center of the box a max imum radial position masr defaults to the smaller value of the box lengths in x and y directions and a minimum and maximum value of z It is possible to also resolve different polar angles thus using it as a fu
154. etes most files that were created during the compilation Will keep for example the built doxygen documentation and the ESPResSo binary dist Creates a tar gz file of the ESPResSo sources This will include all source files as they currently are in the source directory i e it will include local changes This is useful to give your version of ESPResSo to other people The variable extra can be used to specify additional files and directories that are to be included in the archive file Example make dist extra myconfig h internal will create the archive file and include the file myconfig h and the directory internal with all files and subdirectories install Install ESPResSo The variables prefix and exec prefix can be used to specify the installation directories otherwise the defaults defined by the configure script are used prefix sets the prefix where all ESPResSo files are to be installed exec prefix sets the prefix where the executable files are to be installed and is required only when there is an architecture specific directory Example make install prefix usr local will install all files below usr local uninstall Uninstalls ESPResSo i e removes all files that were installed during make install The variables are identical to the variables of the install target ug Creates the User guide in the doc ug subdirectory only when using the develop ment sources dg Creates the Developers guide in the doc dg subdirectory only
155. etting up non equilibrium MD o 6 4 cellsystem Setting up the cell system 6 5 AdResSl 66 CUDAl 6 7 Creating bonds when particles collide 6 8 Catalytic Reactions 6 9 Galilei Transform and Particle Velocity Manipulation 7 Running the simulation 7 1 integrate Running the simulation 7 2 change_volume Changing the box volume ba a a e eo a A E See eee gated ap Gis Ges ee does Gon a Goatees 7 6 Parallel tempering 7 7 Metadynamics 8 Analysis in Tcl 8 1 Available observables 8 2 Analyzing groups of particles molecules o o 8 3 Storing configurations 8 4 uwerr Computing statistical errors in time series Q Analysis in the core 9 1 Observables 9 2 Correlations 10 Input Output 10 1 No generic checkpointing 10 2 blockfile Using the structured file format 10 3 Writing VTF files 10 4 writevtk Particle Visualization in paraview 10 5 Writing PDB PSF files 10 6 Online visualisation with VMD 10 7 Errorhandling 1 Auxilliary commands 11 1 Finding particles and bonds 11 2 Additional Tcl math functions 74 74 76 80 81 82 83 83 85 87 89 89 89 90 90 90 91 95 97 97 107 112 113
156. etween particle 106 and 211 with additional loose particles 212 213 e g counter ions Note that the part command stores any bonds only with the particle of lower particle number which is why part 109 would only return bonds O 110 therefore not revealing the bond between particle 109 and the preceding particle 108 while countBonds would return all bonds particle 109 participates in 11 1 2 findPropPos Syntax findPropPos particle ropertyjist property Description Returns the index of property within particle ropertyjist which is expected to have the same format as part particle d If property is not found 1 is returned This function is useful to access certain properties of particles without hard wiring their index position which might change in future releases of part Example lindex part i findPropPos part i type This returns the particle type id of particle without fixing where exactly that informa tion has to be in the output of part i 139 11 1 3 findBondPos Syntax findBondPos particle roperty ist Description Returns the index of the bonds within particle roperty ist which is expected to have the same format as part particle_number hence its output is the same as findPropPos particle roperty ist bonds If the particle does not have any bonds 1 is returned 11 1 4 timeStamp Syntax timeStamp path prefix postfix suffix Description Modifies the filenam
157. f chains Syntax 1 analyze lt gi gt lt g2 gt lt g3 gt chain_start n_chains chain_length 2 analyze g123 init chain_start n_chains chain_length Description Variant 1 returns e the mean square displacement of the beads in the chain lt g1 gt e the mean square displacement of the beads relative to the center of mass of the chain lt g2 gt e or the motion of the center of mass lt g3 gt averaged over all stored configurations see section 8 3 on the next page At short time scales g1 and g2 coincide since the motion of the center of mass is much slower At large timescales g1 and g3 coincide and correspond to the center of mass motion while g2 levels off g2 and g3 together correspond to g1 For details see Grest and Kremer 21 Variant 2 returns all of these observables for the current configuration as compared to the reference configuration The reference configuration is set when the option init is used Output format variant 1 gi 0 x dt gi 1 dt Output format variant 2 g1 t g2 t g3 t 111 8 3 Storing configurations Some observables i e non static ones require knowledge of the particles positions at more than one or two times Therefore it is possible to store configurations for later analysis Using this mechanism the program is also able to work quasi offline by successively reading in previously saved configurations and storing them to perform any an
158. for the area per lipid to be used in a making a rough calculation of the area of clusters Calls the espresso command analyze aggregation which groups molecules in the sys tem into aggregates Output to time_vs_clust is maximum cluster size minimum cluster size average size of clusters including those of size 2 or greater standard de viation of clusters including those of size 2 or greater number of clusters of size 2 or greater total average cluster size total cluster size standard deviation total number of clusters length of the interface between clusters standard deviation of the interface length number of clusters for which length was calculate Additionally at each call of print_averages the complete size histogram is printed to a file with the formatted name sizehisto format 05d time density_profile nbins arg hrange arg beadtypes arg colloidmoltypes arg r arg nogrid verbose output av_zprof e nbins 100 Number of slices into which the height range is divided for the purpose of calculating densities 163 e hrange 6 The maximum vertical distance from the bilayer midplane for which to calculate densities Note that the complete vertical range is therefore 2 varhrange e beadtypes 0 A tcl list of the bead types for which to calculate a density profile e colloidmoltypes A tcl list of molecule types identifying the molecules which are colloids in the system The default value is a null list e r 0
159. formula has a product decomposition and can be evaluated hierarchi cally similarly to the fast multipole methods For particles sufficiently separated in the z axis one can Fourier transform the potential along both x and y We obtain the far formula as e2 fpaz e2 fpa z z foa FH 1 A P E y 2 Untty D E2TiUy TY e Tiuspe y 2TUzUy ue f p q 0 where Azy z are the box dimensions fpg y Usp uyg fp UP fg Ur Wp 27Ugp and wq 27uyq The advantage of this formula is that it allows for a product decomposition into components of the particles For example e faz e fra 2 2 e2 fpa i 27 fp925 etc Therefore one just has to calculate the sum over all these exponentials on the left side and on the right side and multiply them together which can be done in O N computation time As can be seen easily the convergence of the series is excellent as long as z is sufficiently large By symmetry one can choose the coordinate with the largest distance as z to optimise the convergence Similar to the Lekner sum we need a different formula if all coordinates are small i e for particles close to each other For sufficiently small uyp and u x we obtain the near formula as a cosh 27 fpqz Qriu ae Piia Vue Y a y et p q gt 0 Pa Au Y Ko 2rusppi Kn 27Uxpp 1 cos 27Upx l p gt 0 Que 2 7 R Qty e iy 1 _1 pen 1 uzsr y r 1 uza Ug El E l Gn lge n2 2 lo
160. fter using vs_auto_relate e Virtual sites can not be placed relative to other virtual sites as the order in which the positions of virtual sites are updated is not guaranteed Always relate a virtual site to a non virtual particle placed in the center of mass of the rigid arrangement of particles 40 e Don t forget to declare the particle virtual in addition to calling vs_auto_relate e In case you know the correct quaternions you can also setup a virtual site using part v virtual 1 quat q vs_relative n d where n is the id of the non virtual particle and d is its distance from the virtual site e In a simulation on more than one CPU the effective cell size needs to be larger than the largest distance between a non virtual particle and its associated virtual sites To this aim you need to set the global variable min_global_cut to this largest distance ESPResSo issues a warning when creating a virtual site with vs_auto_relate_to and the cutoff is insufficient e If the virtual sites represent actual particles carrying a mass the inertia tensor of the non virtual particle in the center of mass needs to be adapted 4 4 3 Additional features The behaviour of virtual sites can be fine tuned with the following switches in myconfig h sec 3 4 e VIRTUAL_SITES_NO_VELOCITY specifies that the velocity of virtual sites is not com puted e VIRTUAL_SITES_THERMOSTAT specifies that the Langevin thermostat should also act on virtua
161. g 47 Note that this time we calculate instead of i e we omit the contribution of the primary simulation box This is very convenient as it includes the case of self energy and makes a smooth function To obtain one has to add the 1 r contribution of the primary box The self energy is given by 1 fq 7fpa z 1 Both the near and far formula are derived using the same convergence factor approach and consequently the same singularity in 6 is obtained This is important since otherwise the charge neutrality argument does not hold To obtain the O N log N scaling some algorithm tricks are needed which are not used in MMM1D MMM2D or ELC and are therefore not discussed here For details see Strebel 51 MMM is not implemented in ESPResSo 8uz y Ky QruzAypl 2u 40 1 2 log 47 l p gt 0 3 0 0 0 2uxtty Y p q gt 0 211 E 2 MMM2D In the case of periodicity only in the x and y directions the far formula looks like 27 fpglz x y z MgUy Y pg gt o COS Wpx cos wyY 2r fq z 2r fpl z 2UgUy sd Az 7 cos wqy yu T cos up r _ 274 Uy 2 and the near formula is O x Yy 2 Aug i p gt 0 Kolwppr Ko wpp 1 cos wpx e AS Ny 1 2s Donor ray ati Dea A 1 gt 5 Ye Ny urxr y2 Ny u1 2 thy nl 3 Y Y n 2n uxp 2uy log ar y As said before the energy obtained from these potentials is equal to the electrostatic en
162. g a checkpoint integrate set_checkpoint integrate has only one call to the force calculation routine while read_checkpoint integrate has two at the beginning of the integrate command because loading a new system from disk typically requires re initializing the system and since the forces routine also uses the thermostat which in turn draws random numbers the two situations do not end up at the same segment of the random number sequence all random events will therefore slightly differ To prevent this simply include a call to invalidate_system upon setting the checkpoint because in that case both scenarios will call the forces routine twice at the beginning of the second integration phase thus having their random number sequences in total sync Without applying this command directly before or after writing a checkpoint you will run into a different state of the random number generator when reading the checkpoint to start again later 90 7 6 Parallel tempering Syntax parallel_tempering main rounds N swap swap perform perform init init values T connect master port port load jnode resrate Nresei info info Description This command can be used to run a parallel tempering simulation Since the simulation routines and the calculation of the swap probabilities are provided by the user the method is not limited to sampling in the temperature space However we assume in the following that the sampled value
163. ge and you have a higher probability that it is answered soon e Your question and the answers are archived and the archives can be searched by others e The answer may be useful also to other registered users e There may not be a unique answer to your problem and it may be useful to get suggestions from different people Please remember that this is a community mailing list It is other users and developers who are answering your questions They do it in their free time and are not paid for doing it 184 16 2 Contributing your own code If you are planning to make an extension to ESPResSo or already have a piece of your own code which could be useful to others you are very welcome to contribute it to the community Before you start making any changes to the code you should obtain the current development version of it For more information about how to obtain the development version refer to the homepage It is also generally a good idea to contact the mailing lists before you start major coding projects It might be that someone else is already working on the problem or has a solution at hand 16 3 Developers guide Besides the User guide ESPResSo also contains a Developers guide which is a program mer documentation automatically built from comments in the source code and using Doxygen It provides a cross referenced documentation of all functions and data struc tures available in ESPResSo source code It can be built
164. guments in their correct order for a particular call to the espresso inter command See the espresso inter command for a list of possible non bonded interactions and correct syntax Set all the bonded interactions mbtools utils init_random n_procs e n_procs The number of processors used in this job Initialize the random number generators on each processor based on the current time with a fixed increment to the time seed used for each proc mbtools utils initialize_vmd flag outputdir ident extracommands arg e flag Depending on the value of this parameter initialize vmd to one of its possible states interactive VMD is started and a connection to espresso established for immediate viewing of the current espresso process With some luck this might even work sometimes If VMD doesn t get a proper connection to espresso then it will crash offline Just constructs the appropriate psf and vmd_animation script files and writes them to the output directory so that pdb files generated with writepdb can be viewed with vmd e vmd_animation script default Any value other than those above for flag will just result in vmd not being initialized e outputdir The directory where vmd output will be written e ident A basename to be be given to vmd files e extracommands A list of strings each of which will be written to the end of the vmd_animationscript Use this to give additional commands to vmd Pr
165. h step a n ko 02 Teut Required features SMOOTH_STEP Description This defines a smooth step interaction between particles of the types typel1 and type2 for which the potential is V r 01 d 1 exp 2ko r 02 5 5 for r lt reut and V r 0 elsewhere With n around 10 the first term creates a short range repulsion similar to the Lennard Jones potential while the second term provides a much softer repulsion This potential therefore introduces two length scales the range of the first term 01 and the range of the second one 02 where in general 01 lt 02 5 1 6 BMHTF potential Syntax inter typel type bmhtf nacl A B C Do Tous Required features BMHTF_NACL Description This defines an interaction with the short ranged part of the Born Meyer Huggins Tosi Fumi potential between particles of the types type and type2 which is often used to 46 simulate NaCl crystals The potential is defined by V r Aexp B o r Cr Dro esits 5 6 where esnir is chosen such that V reut 0 For r gt few the V r 0 For NaCl the parameters should be chosen as follows types A kJ mol B A 1 C ASkJ mol D A8kI mol oO A Na Na 25 4435 3 1546 101 1719 48 1771 2 34 Na Cl 20 3548 3 1546 674 4793 837 0770 2 755 Cl1 Cl 15 2661 3 1546 6985 6786 14031 5785 3 170 The cutoff can be chosen relatively freely because the potential decays fast a value around 10 seems reasonab
166. h value lindex rdf 1 lappend rlist lindex value 0 lappend rdflist lindex value 1 The shown analyze rdf command returns the distribution function of particles of type 1 around particles of type 0 i e of opposite charges for radii between 0 9 and half the 16 Figure 2 1 VMD Snapshot of the salt system box length subdivided into 100 bins Changing the first two parameters to either 0 0 or 1 1 allows to determine the distribution for equal charges The result is a list of r and g r pairs which the following foreach loop divides up onto two lists rlist and rdflist To average over a set of configurations we put the two last code snippets into a loop like this set cnt 0 for set i 0 i lt 100 incr i lappend avg_rdf 0 foreach filename argv set f open filename r while blockfile f read auto eof close f set rdf analyze rdf 0 1 0 9 expr box_1 2 100 set rlist set rdflist foreach value lindex rdf 1 lappend rlist lindex value 0 lappend rdflist lindex value 1 set avg_rdf vecadd avg_rdf rdflist incr cnt set avg_rdf vecscale expr 1 0 cnt avg_rdf Initially the sum of all g r which is stored in avg_rdf is set to 0 Then the loops over all configurations given by argv calculates g r for each configuration and adds up all the g r in avg_rdf Finally this sum is normalized by dividing by the number of configurations Note the
167. he averaged internal distances within the chains over all pairs of particles If lt internal_dist gt is used the values are averaged over all stored configurations see section 8 3 on page 112 Output format idf 0 idf 1 idf chain_length 1 y The index corresponds to the number of beads between the two monomers considered 0 next neighbours 1 one monomer in between Internal distances II specific monomer Syntax analyze bond_dist lt bond_dist gt index inder chain_start n_chains chain_length Description In contrast to analyze internal_dist it does not average over the whole chain but rather takes the chain monomer at position index default 0 i e the first monomer on the chain to be the reference point to which all internal distances are calculated If lt bond_dist gt is used the values will be averaged over all stored configurations see section 8 3 on page 112 Output format bdf 0 bdf 1 bdf chain_length 1 index 109 Bond lengths Syntax analyze bond_1 lt bond_1 gt chain_start n_chains chain_length Description Analyzes the bond lengths of the chains in the system Returns its average the standard deviation the maximum and the minimum If you want to look only at specific chains use the optional arguments i e chain_start 2x MPC and n_chains 1 to only include the third chain s monomers If lt bond_1 gt is used the value will be
168. he connectivity of the particle up to range defaults to 1 For particle 5 in a linear chain the result up to range 3 would look like 28 4 6 L433 67 432 678 The function is useful when you want to create bonded interactions to all other particles a certain particle is connected to Note that this output can not be used as input to the part command Check results if you use them in ring structures If none of the options is specified it returns all properties of the particle if it exists in the form O pos 2 1 6 4 3 1 type 0 q 1 0 v 0 0 0 0 0 0 f 0 0 0 0 0 0 bonds 0 480 0 368 which may be used as an input to this function later on The first integer is the particle number Variant 2 returns the properties of all stored particles in a tcl list with the same format as specified above 0 pos 2 1 6 4 3 1 type 0 q 1 0 v 0 0 0 0 0 0 f 0 0 0 0 0 0 bonds 0 480 0 368 1 pos 1 0 2 0 3 0 type 0 q 1 0 v 0 0 0 0 0 0 f 0 0 0 0 0 0 bonds 0 340 0 83 2 3 4 1 3 Deleting particles Syntax 1 part pid delete 2 part deleteall Description In variant 1 the particle pid is deleted and all bonds referencing it Variant 2 will delete all particles currently present in the simulation Variant 3 will delete all currently defined exclusions 4 1 4 Exclusions Syntax 1 part auto_exclusions range 2 part delete_exclusions Required features
169. he object immersed in the fluid will move under the influence of the deforming forces defined through the bonds and under the influence of the fluid motion This interaction is based on the frictional force between the fluid and the surface particles Therefore the object will move in the flow only if there is a nonzero difference between the fluid velocity and the particle velocity In other words there has to be at least small flow through the membrane which is in most cases unphysical However this unphysical flow through the membrane is probably negligible in larger scales 154 13 1 Membranes With this approach it is easy to model also elastic sheets or free membranes that do not necessarily enclose a 3D object In this case area_force_global and volume_force interactions are not needed since these two interactions are ment for closed immersed objects 13 2 Parameters There are several parameters involved in this model All of them should be calibrated according the application e Mass of the particles Every particle has its mass that influences the dynamics e Friction coefficient The main parameter describing the fluid particle interaction is the friction parameter from the ESPResSo command 1bfluid e Parameters of elastic moduli Elastic behaviour can be described by five different eleastic moduli stretching bending local and global area preservation and volume preservation Each of them has its own scaling parame
170. her or not one or some of the features are compiled into the current program with help of the following Tcl commands e code_info provides information on the version compilation status and the debug status of the used code It is highly recommended to store this information with your simulation data in order to maintain the reproducibility of your results Exemplaric output ESPRESSO vi 5 Beta Neelix Last Change 23 01 2004 Compilation status PARTIAL_PERIODIC ELECTROSTATICS EXTERNAL_FORCES CONSTRAINTS TABULATED LENNARD_JONES BOND_ANGLE_COSINE Debug status MPI_CORE FORCE_CORE e has_feature lt feature gt 146 tests if feature is compiled into the ESPResSo kernel A list of possible features and their names can be found here require_feature lt feature gt tests if feature is feature is compiled into the ESPResSo kernel will exit the script if it isn t and return the error code 42 A list of possible features and their names can be found here 147 12 Lattice Boltzmann For an implicit treatment of a solvent ESPResSo allows to couple the molecular dynam ics simulation to a Lattice Boltzmann fluid The Lattice Boltzmann Method LBM is a fast lattice based method that in its pure form allows to calculate fluid flow in different boundary conditions of arbitrarily complex geometries Coupled to molecular dynamics it allows for the computationally efficient
171. his LUB requirements or will throw an error if there are none The gap size gap_size gives the height of the empty region between the system box and the neighboring artificial images again see the paper ESPResSo does not make sure that the gap is actually empty this is the users responsibility The method will compute fine of the condition is not fulfilled however the error bound will not be reached Therefore you should really make sure that the gap region is empty e g by constraints The setting of the far cutoff far_cutoff is only intended for testing and allows to directly set the cutoff In this case the maximal pairwise error is ignored The periodicity has 66 to be set to 1 1 1 still and the 3d method has to be set to epsilon metallic i e metallic boundary conditions For details see appendix E on page 210 By default ELC just as P3M adds a homogeneous neutralizing background to the system in case of a net charge However unlike in three dimensions this background adds a parabolic potential across the slab 8 Therefore under normal circumstance you will probably want to disable the neutralization using noneutralization This corresponds then to a formal regularization of the forces and energies 8 Also if you add neutralizing walls explicitely as constraints you have to disable the neutralization The dielectric contrast features work exactly the same as for MMM2D see the docu mentation above Make sure
172. hree body potentials are implemented following the procedure in Ref 53 A different formula is used to calculate contribution from electrostatic interac tions in P3M For electrostatic interactions the k space contribution is not well tested so use with caution Anything outside that is currently not implemented Four body dihedral potentials are not included In case of rigid body rotation virial contribution from torques is not included Constraints of any kind are not currently accounted for in the pressure calculations The pressure is no longer correct e g when particles are confined to a plane The command is implemented in parallel Output format variant 1 pressure total_pressure ideal ideal_gas pressure 4 bond_type pressure 4 nonbonded_type pressure coulomb pressure specifying the pressure the ideal gas pressure the contributions from bonded interac tions the contributions from non bonded interactions and the electrostatic contributions 8 1 22 Stress Tensor Syntax 1 analyze stress_tensor 2 analyze stress_tensor total 3 analyze stress_tensor totals ideal coulomb tot_nonbonded_inter tot_nonbonded_intra 4 analyze stress_tensor bonded bond ype 5 analyze stress_tensor nonbonded typeidl typeid2 6 analyze stress_tensor nonbonded_intra typeid 7 analyze stress_tensor nonbonded_inter typeid 105 Description Computes the stress tensor of the syste
173. i gt jels S lt S 1 qje T tea cos wpztj sin wqyj go sin wpxi cos wqyi J jeIs S lt S 1 qje Tra i sin Wpx cos wgyj a sin Wpx sin WyYi gt jers s lt s 1 qje tra sin Wpx sin wyy 212 This representation has the advantage that the contributions of the two particles are decoupled For all particles j only the eight terms ei qye 2 624 sin cos wpr sin cos wyy are needed The upper index describes the sign of the exponential term and whether sine or cosine is used for x and yj in the obvious way These terms can be used for all expressions on the right hand side of the product decomposition Moreover it is easy to see from the addition theorem for the sine function that these terms also can be used to calculate the force information up to simple prefactors that depend only on p and q E s c s c and adds them up Every processor starts with the calculation of the terms 33 in each layer so that one obtains s c s c _ 4 8 8 c s u jJESs 1 Now we calculate El 8 e s c y a 05 0 t lt s 1 and h s c s c 8 c 5 c an gt t gt s l which are needed for the evaluation of the product decomposition While the bottom l s c s c processor can calculate s directly the other processors are dependent on its results Therefore the bottom processor starts with the calculation of its alls 8 and sends up bs 8 and o o of
174. ible to build a single binary that can satisfy all needs A user should always activate only those features that are actually needed This means however that learning how to compile ESPResSo is a necessary evil The build system of ESPResSo uses the GNU autotools which are developed since more than 20 years and allow to compile software easily on a wide range of platforms 3 1 Running configure The first step of building ESPResSo is to run the shell script configure which is to be found in the top level source directory The script collects all the information required by the compilation process It will determine how to use and where to find the compiler as well as the different libraries and tools required by the compilation process and it will test what compiler flags are to be used The script will find out about most of these things automatically If something is missing it will complain and give hints how to solve the problem The generic syntax of calling the configure script is configure options variable value http espressomd org http git org https savannah nongnu org projects espressomd 20 If you are using the development source code from the git repository before you can call configure it is necessary to have the GNU autotools autoconf and automake installed Then you can call the script bootstrap sh from the top level source directory which will generate the configure script 3
175. ic field in x direction by half a time step 13 a pdate the particle positions in x direction by half a time step E 14 Update the B field by half a time step a 15 Update the particle momenta by half a time step D 5 Self energy The interpolation of the charges onto the lattice gives rise to the artificial force exerted on the particle by its own field In order to cure this remedy the direct subtraction of the self energy is introduced For the interpolated charge cloud the self energy can be directly calculated For the simple cubic lattice in three dimensions the linear interpolation will give 8 charges which are placed at the corners of the cube with edge length a see Fig D 3 q dg Pi q UE ab q N a 9 Figure D 3 Linear interpolation scheme 207 Therefore in our case the self energy is a symmetric bilinear form defined by the matrix a the elements of which do not depend on the position of the charge In our algorithm the values of the coefficients are 3 cca R B 3 D11 Qij e d aan 53 1 cos k where L is the number of lattice points per dimension R coordinates of the interpo lated charges and k the wave vector Those values are calculated during the initialization step and are used in the calculation of the self force The value of the self force which has to be subtracted from the overall forces is given by the following ansatz mt aoe o 2
176. icle and bonding information t arg The topology file corresponding to the file to be read tol arg 0 000001 Tolerance for comparison of box dimensions Use particle positions contained in a file to initialise the locations of particles for a particular geometry The box dimensions in the file and those set by the user are compared and an error is returned if they are not the same to within a tolerance value of tol Even though we read from a file we also generate a topology from the ny olslist and this topology is compared with the topology that is read in to check if the number of particles are the same singlemol c arg o arg trapflag arg ctrap arg trapspring arg bondl arg c arg 0 00 0 0 0 The molecule center Exactly what this means depends on the molecule type e o arg 0 00 0 1 0 The orientation vector for the molecule This is also molecule type dependent e trapflag arg 000 Set this optional argument to cause a molecule to be trapped by its center of mass You should give three integers corresponding to each of the three coordinate axes If a value of 1 is given then motion in that axis is trapped e ctrap arg Set this optional argument to the central point of the trap This works much like an optical trap in that molecules will be attracted to this point via a simple harmonic spring force e trapspring arg 20 The spring constant for the trap potential harmonic spring Simply place a singl
177. icles This also applies to the electrostatic interaction however due to its long ranged nature it requires special care and ESPResSo handles it separately with a number of state of the art algorithms The particle type and the charge are both defined using the part command A bonded interaction defines an interaction between a number of specific particles it only applies to the set of particles for which it has been explicitely set A bonded interaction between a set of particles has to be specified explicitely by the part bond command while the inter command is used to define the interaction parameters Syntax inter Description Without any arguments inter returns a list of all defined interactions as a Tcl list The format of each entry corresponds to the syntax for defining the interaction as described below Typically this list looks like 0 0 lennard jones 1 0 2 0 1 1225 0 0 0 0 0 FENE 7 0 2 0 5 1 Isotropic non bonded interactions Syntax inter typel type2 interaction parameters Description This command defines an interaction of type interaction between all particles of type typel and type2 The possible interaction types and their parameters are listed below If the interaction is omitted the command returns the currently defined interaction between the two types using the syntax to define the interaction e g O O lennard jones 1 0 2 0 1 1225 0 0 0 0 For many non bonded interactions it is possible to
178. icles in the system if no type is given Returns a Tcl list containing the squared radius of gyration three shape descriptors asphericity acylindricity and relative shape anisotropy eigen values of the gyration tensor and their corresponding eigenvectors The eigenvalues are sorted in descending order 8 1 17 Aggregation Syntax analyze aggregation dist_criteria s_mol_id f_mol_id min_contact charge_criteria Description Returns the aggregate size distribution for the molecules in the molecule id range s_mol_id to f_mol_id If any monomers in two different molecules are closer than dist_criteria they are considered to be in the same aggregate One can use the op tional min_contact parameter to specify a minimum number of contacts such that only molecules having at least min_contact contacts will be considered to be in the same aggregate The second optional parameter charge_criteria enables one to consider ag gregation state of only oppositely charged particles 8 1 18 Identifying pearl necklace structures Syntax analyze necklace pearl_threshold back_dist space_dist first length Description Algorithm for identifying pearl necklace structures for polyelectrolytes in poor solvent The first three parameters are tuning parameters for the algorithm pearl_threshold is the minimal number of monomers in a pearl back_dist is the number of monomers along the chain backbone which are excluded from the space distance criterion to fo
179. id which acts between two particles and actually subtracts the Lennard Jones interaction between the involved particles The first parameter reserved is a dummy just kept for compatibility reasons The second parameter R is used as a check if any bond length in the system exceeds this value the program terminates When using this interaction it is worthwhile to consider capping the Lennard Jones potential appropriately so that round off errors can be avoided This interaction is useful when using other bond potentials which already include the short ranged repulsion This often the case for force fields or in general tabulated potentials 5 3 4 Rigid bonds Syntax inter bondid rigid_bond constrained_bond_distance positional_tolerance velocity_tolerance Description To simulate rigid bonds ESPResSo uses the Rattle Shake algorithm which satisfies inter nal constraints for molecular models with internal constraints using Lagrange multipliers 53 5 3 5 Tabulated bond interactions Syntax 1 inter bondid tabulated bond filename 2 inter bondid tabulated angle filename 3 inter bondid tabulated dihedral filename Description This creates a bond type with identificator bondid with a two body bond length variant 1 three body angle variant 2 or four body dihedral variant 3 tabulated poten tial The tabulated forces and energies have to be provided in a file filename which is formatted identically as the files
180. iding the details of the complex numerical algorithms The Tcl interpreter contains several special commands as an extension to Tcl which provide the interaction with the simulation engine Thus the user has at hand the full realm of Tcl commands and constructs plus a few commands to communicate with the simulation engine The interfacing commands are designed so that they can both set properties of the system set up particles interactions thermostat and retrieve information about the already set up entities The standard Tcl constructs allow for a flexible decision making in the course of the simulation This can be for example exploited to check whether a simulation has reached the desired state With a certain overhead in efficiency it can also be used to reject accept new configurations in combined MD MC schemes In principle any parameter which is accessible from the Tcl level can be changed at any moment of runtime In this way methods like thermodynamic integration become readily accessible The focus of the user guide is documenting the Tcl commands their behaviour and use in the simulation It only describes certain technical details of implementation which are necessary for understanding how the commands work Technical documentation of the code and program structure is contained in the Developers guide see section 16 3 1 4 On units What is probably one of the most confusing subjects for beginners of ESPResSo is that ESPR
181. ij and gt Mij There are two implemented cases DPD_MASS_RED uses the reduced mass MiMi Mi y pi ea Mi Mj while DPD_MASS_LIN uses the real mass mass Mi Mj Mij 2 The prefactors are such that equal masses result in a factor 1 79 6 2 4 Isotropic NPT thermostat Syntax thermostat npt_isotropic temperature gamma0 gammaV Required features NPT Description This theormstat is based on the Anderson thermostat see 1 34 and will thermalize the box geometry It will only do isotropic changes of the box Be aware that this feature is neither properly examined for all systems nor is it main tained regularly If you use it and notice strange behaviour please contribute to solving the problem 6 3 nemd Setting up non equilibrium MD Syntax 1 2 nemd exchange n_slabs n_exchange 2 nemd shearrate n_slabs shearrate 3 nemd off 4 nemd 5 6 nemd viscosity nemd profile Required features NEMD Description Use NEMD Non Equilibrium Molecular Dynamics to simulate a system under shear with help of an unphysical momentum change in two slabs in the system Variants 1 and 2 will initialise NEMD Two distinct methods exist Both methods divide the simulation box into n_slab slabs that lie parallel to the x y plane and apply a shear in x direction The shear is applied in the top and the middle slabs Note that the methods should be used with a DPD thermostat or in an NVE ensemble Fu
182. in the simulation by setting a list called moltypes In this case two different lipid types are setup and assigned to moltypeids 0 and 1 respectively Moltype 0 will consist of three beads per lipid the first of which is of atomtype 0 and the second and third of which are of atomtype 1 Bonds in the lipid will be of type 0 and 1 respectively see the mbtools system_ generation place_lipid_linear function for further details set moltypes list O lipid 011 01 i lipid 0222 4 02 We then construct system specs for a flat bilayer and a spherical bilayer and group these into a system_specs list First the spherical system_specs 167 set geometry geometry sphere shuffle c 0 0 0 0 15 0 set n_molslist n_molslist 0 1000 lappend spherespec geometry lappend spherespec n_molslist The flat system_spec set geometry geometry flat fixz set n_molslist n_molslist 1 3000 lappend bilayerspec geometry lappend bilayerspec n_molslist Now group together the system pecs into a master list lappend system_specs spherespec lappend system_specs bilayerspec Make the call to setup_system mbtools system_generation setup_system system_specs setmd box_1 moltypes mbtools system_generation get_trappedmols returns the internal list variable trappedmols which keeps track of all molecules that have been trapped by their center of mass This function should be called af
183. ined by its z midplane zo its total thickness 6z and the interface width Kk Therefore the interaction strength is e everywhere except for the region of the box zo 02 2 lt z lt zo 02 2 The interface width smoothly interpolates between the two regions to avoid discontinuities As an example one can think of modeling hydrogen bonds in two different environments water where the interaction is rather weak and in a lipid bilayer where it is comparatively stronger 5 2 2 Gay Berne interaction Syntax inter typel type gay berne 9 Oo Teutor kl k2 wv Required features ROTATION GAY_BERNE Description This defines a Gay Berne potential for prolate and oblate particles between particles of the types type1 and type2 The Gay Berne potential is an anisotropic version of the classic Lennard Jones potential with orientational dependence of the range og and the well depth g Assume two particles with orientations given by the unit vectors and and inter molecular vector r rf If r lt rey then the interaction between these two particles is given by V rij aj s de Fiz Gj s Ti a 5 5 15 otherwise V r 0 The reduced radius is r o f Uy 00 ki 5 16 00 2 27173 1 FU Pa t o f 00 1 X u ue a A u Uy 5 17 2 1 x j 1 xu j i Og cee ey ef a ee eee eT ey Ee eo 1 x a aj 1 4 cm oe y T 2 4 5 18 2 1 a a 1
184. ion Copy a group of particles including their bonds Positions can be shiftled by an offset s_x s_y s_z otherwise the copied set is at exactly the same position as the original set The particles can be given as a combination of list s or range s The new particles obtain in any case consecutive identities after the largest current identity Bonds within the defined particle set are copied with translated identities but not bonds with particles 35 outside the list That is if the particle set corresponds to a molecule intramolecular bonds are preserved but not intermolecular ones Examples of use copy_particles set 1 2 3 4 shift 0 0 0 0 0 0 copy_particles set 1 2 set 3 4 copy_particles range 1 4 All these examples do the same making exact copies of particles 1 through 4 4 3 constraint Setting up constraints Syntax 1 constraint wall normal nz ny n dist d type id penetrable flag reflecting flag 2 constraint sphere center Cy cy Cz radius rad direction direction type id penetrable flag reflecting flag 3 constraint cylinder center Cyr Cy Cz axis Nz Ny Nn radius rad length length direction direction type id penetrable flag reflecting flag 4 constraint rhomboid corner Pz Py Pz a dy dy Gz b by by bz C Cy Cy Cz direction direction type id penetrable flag reflecting flag 5 constraint maze nsphere n dim d sphrad r cylrad r type id penetrable flag 6 constraint pore center Cy Cy Cz axis Nz N
185. ions verbose output powav dat Routine for calculating the power spectrum of height and thickness fluctuations for a flat bilayer sheet Uses the modes_2d routine in ESPResSo to calculate the height and thick ness functions and perform the fft See the documentation in the file fluctuations tcl for detail on what is calculated and how to obtain a stiffness value from the resulting output Note that this routine causes a crash if it detects a large hole in the bilayer localheights range arg nbins arg rcatch arg verbose output av_localh e range 1 0 Range of local height deviations over which to bin 164 e nbins 100 Number of slices to divide up the height range into for the purposes of creating a profile e rcatch 1 9 The distance about a single lipid to use a starting value for finding the 6 closest neighbours For each lipid we calculate its 6 nearest neighbours and then calculate the height differ ence between the central lipid and these neighbours Taking these 6 values for each lipid we then create a histogram of number densities as a function of the height difference localorients range arg nbins arg verbose output av_localo e range 1 0 Range of orientation deviations to consider e nbins 100 Number of bins to use for histogram Calculates the projection of the lipid orientation vector onto the xy plane for each lipid and then bins the absolute values of these vectors orient_order ver
186. iour with respect to statistical and systematic errors please read the cited literature This type of correlator has been in use for years in the analysis of dynamic light scattering 45 About a decade later it found its way to the Fluorescence Correlation Spectroscopy FCS 32 The book of Frenkel and Smit describes its application for the special case of the velocity autocorrelation function Let us consider a set of N observable values as schematically shown in Figures where a value of index i was measured in time idt We are interested in computing the correlation function according to Equation for a range lag times 7 i 7 dt between the measurements and j To simplify the notation we further drop t when referring to observables and lag times The trivial implementation takes all possible pairs of values corresponding to lag times T Tmin Tmax Without loss of generality let us further consider Tmin 0 The computational effort for such an algorithm scales as O T ax As a rule of thumb this is feasible if Tmax lt 10 The multiple tau correlator provides a solution to compute the correlation functions for arbitrary range of the lag times by coarse graining the high 7 values It applies the naive algorithm to a relatively small range of lag times T 0 p 1 This we refer to as compression level 0 To compute the correlations for lag times 7 p 2 p 1 the original data are first coarse grained so that
187. iption Returns a Tcl list of the particle ids of all particles within a given radius r_catch around the position of the particle with number pid in variant 1 or around the spatial coordi nate 1 y z in variant 2 8 1 3 Particle distribution Syntax analyze distribution part_type_list_a part_type_list_b rmin rmaz rbins log_flag int_flag Description Returns its parameters and the distance distribution of particles with types specified in part_type_list_a around particles with types specified in part_type_list_b with distances between rmin and rmaz binned into rbins bins The bins are either equidistant if log_flag 0 or logarithmically equidistant if log_flag gt 1 If an integrated distribution is required use int_flag 1 The distance is defined as the minimal distance between a particle of one group to any of the other group Output format The output corresponds to the blockfile format see section 10 2 on page 127 parameters r dist r 8 1 4 Radial density map Syntax analyze radial_density_map zbins ybins xrange yrange axisofrotation centerofrotation beadtypelist thetabins Description Returns the radial density of particles around a given axis Parameters are e rbins histogram bins in x direction e ybins histogram bins in y direction 98 e rrange range for analysis in x direction e yrange range for analysis in y direction e azxisofrotation rotate around given axi
188. is set to on gamma_0 Required if npt is on Usually set to 1 as for langevin gamma gamma_v Required if npt is on Box friction use_umd offline vmd mode mgrid 8 The number of meshpoints per side for dividing the bilayer plane into a grid stray_cut_off 1000 0 Distance of the end tail bead from the bilayer midplane beyond which a lipid is deemed to have strayed from the membrane bulk systemtemp The temperature of the simulation during the main run outputdir Directory for output tabledir Directory where forcetables are kept ident a name for the simulation topofile the name of the file where the topology is written Usually ident top tablenames A list of forcetable names to be used setbox_l Box dimensions bonded_parms A complete list of the bonded interactions required nb_interactions A complete list of the non bonded interactions required system_specs A list of system specifications see documentation for the setup_ system command in 14 5 moltypes A list of molecule types see documentation in 14 5 warm_time_step timestep to be used during warmup integration main_time_step timestep for the main integration run verlet_skin skin used for verlet nesting list criterion langevin_gamma langevin friction term int_n_times number of times to do main integration int_steps number of steps in each main integration 161 e analysis_write_frequency How often to calculate the analysis e write_frequency
189. izes This method is called electrostatic layer correction ELC The advantage of this approach is that for the image layers z is necessarily large enough so that all interactions can be represented using the product decomposition This allows for an order N evaluation of the ELC term The electrostatic layer correction term is given by N Ete Y gt qigyh pi py ij 1 where h 27 x y z AugUy ae O cos w x cos wyy cosh 27 fpz 24 Uy are Felz h 27 f 242 Uy go A cos Wwyy cos w 2 The implementation is very similar to MMM2d except that the separation between slices closeby and above and below is not necessary E 5 Errors Common to all algorithms of the MMM family is that accuracy is cheap with respect to computation time More precisely the maximal pairwise error i e the maximal error of the y expression decreases exponentially with the cutoffs In turn the computation time grows logarithmically with the accuracy This is quite in contrast to the Ewald methods for which decreasing the error bound can lead to excessive computation time For example P3M cannot reach precisions above 1075 in general The precise form of the error estimates is of little importance here for details see Arnold et al 6 216 One important aspect is that the error estimates are also exponential in the non periodic coordinate Since the number of closeby and far away particles is different for particles near the borde
190. ject These two files are sufficient to describe the geometry and topology of the triangu lation For the definition of bonded interactions the following geometric entities are necessary position of the particles edges lengths of the edges triangles areas of tri angles angles between two triangles sharing a common edge surface of the immersed object volume of the immersed object All these geometrical entities can be computed using the information from the files mesh nodes dat and mesh triangles dat and the computation is done in the script scripts object_in_fluid tcl The script scripts object_in_fluid tcl reads both mesh files generates list of edges and computes all geometrical entities needed for definition of bonded interactions It then executes commands creating the particles interactions and bonds An example of part command is as follows part 0 pos 3 0 3 0 6 0 type 1 mol 1 mass 1 Note that there is feature mol used for the particles With this feature we distinguish between different objects The limit for the number of objects is 10000 However it can be increased by changing the MAX_OBJECTS_IN_FLUID constant Next example shows an interaction inter 106 stretching_force 4 6 5 0 156 By this command btw this command is invisible for the user it is execute inside the scripts object_in_fluid tcl script an interaction for stretching is created with ID 106 Detailed description of the available types of intera
191. l sites e THERMOSTAT_IGNORE_NON_VIRTUAL specifies that the thermostat does not act on non virtual particles 4 5 Grand canonical feature For using ESPResSo conveniently for simulations in the grand canonical ensemble or other purposes when particles of certain types are created and deleted frequently Syntax part gc type find delete status number type Required features GRANDCANONICAL Description By giving only a particle type as an argument ESPResSo will initialize arrays on the c level which will keep track of particles of the given type When using the keyword find and a particle type the command will return a randomly chosen particle id for a particle of the given type Similarly using the delete variant will delete a randomly chosen particle of the given type Notice however that changing the type of an existing particle will not change the arrays If the change is permanent then re initializing the 41 arrays of the involved types will remove this inconsistency The keyword status will return a list with all particles with the given type similarly giving number as argument will return the number of particles which sharing the given type 42 5 Setting up interactions In ESPResSo interactions are set up and investigated by the inter command There are mainly two types of interactions non bonded and bonded interactions Non bonded interactions only depend on the type of the two involved part
192. latter one as temperature However the simulation will crash ESPResSo complains about particle coordinates being out of range The reason for this is simple Due to the initial random setup the overlap energy is around a million kT which we first have to remove from the system In ESPResSo this is can be accelerated by capping the forces i e modifying the Lennard Jones force such that it is constant below a certain distance Before the integration loop we therefore insert this equilibration loop for set cap 20 cap lt 200 incr cap 20 puts t setmd time E analyze energy total inter forcecap cap integrate integ_steps inter forcecap 0 This loop integrates the system with a force cap of initially 20 and finally 200 The last command switches the force cap off again With this equilibration the simulation script runs fine However it takes some time to simulate the system and one will probably like to write out simulation data to configuration files for later analysis For this purpose ESPResSo has commands to write simulation data to a Tcl stream in an easily parsable form We add the following lines at end of integration loop to write the configuration files config_0 through config_19 There also exists a Tcl function degrees_of_freedom which does the same 15 set f open config i w blockfile f write tclvariable box_l density blockfile f write variable box_1 blockfile f wri
193. le In addition to this short ranged interaction one needs to add a Coulombic long ranged part If one uses elementary charges i e a charge of q 1 for the Na particles and q 1 for the Cl particles the corresponding prefactor of the Coulomb interaction is 1389 3549A kJ mol 5 1 7 Morse interaction Syntax inter type type morse Q Tmin Teut Required features MORSE Description This defines an interaction using the Morse potential between particles of the types type and type2 It serves similar purposes as the Lennard Jones potential but has a deeper minimum around which it is harmonic This models the potential energy in a diatomic molecule This potential allows capping the force using inter forcecap see section 5 9 5 For r lt Teut this potential is given by V r e exp 2a r rmin 2exp a r rmin esnitt 5 7 where snift is again chosen such that V reut 0 For r gt Tens the V r 0 5 1 8 Buckingham interaction Syntax inter type type buckingham A B C D Teut Tdiscont shift Required features BUCKINGHAM Description This defines a Buckingham interaction between particles of the types type and type2 for which the potential is given by V r Aexp Br Or Dir eshift 5 8 47 for Tdiscont lt Y lt Teut Below raiscont the potential is linearly continued towards r 0 similarly to force capping see below Above r reyt the
194. le Lennard Jones particle liquid Shows the basic features of ESPResSo How to set up system parameters particles and interactions How to warm up and integrate How to write parameters configurations and observables to files How to handle the connection to VMD pe_solution tcl Polyelectrolyte solution under poor solvent condition Test case for com parison with data produced by polysim9 from M Deserno Note that the equili bration of this system takes roughly 150007 pe_analyze tcl Example for doing the analysis after the actual simulation run offline analysis Calculates the integrated ion distribution P r for several different time slaps compares them and presents the final result using gnuplot to generate some ps files harmonic_oscillator tcl A chain of harmonic oscillators This is a T 0 simulation to test the energy conservation espresso_logo tcl The ESPResSo logo the exploding espresso cup has been created with this script It is a regular simulation of a polyelectrolyte solution It makes use of some nice features of the part command see section 4 1 on page 26 namely the capability to fix a particle in space and to apply an external force 203 D Maxwell Equations Molecular Dynamics MEMD In this chapter we want to give a more thorough introduction to the MEMD or Maggs algorithm for the calculation of Coulomb interactions that is implemented in ESPResSo For an even more detailed description we r
195. le and some data The data itself may consist again of such blocks 127 An example is file Demonstration of the block format variable epsilon _dval_ 1 gt variable p3m_mesh_offset _dval_ 5 0000000000e 01 5 0000000000e 01 5 0000000000e 01 variable node_grid _ival_ 2 2 2 end Whitespace will be ignored within the format space tab and return The keyword variable should be used to indicate that a variable definition follows in the form name data data itself is a block with title _ival_ or _dval_ denoting integer resp double values which then follow in a whitespace separated list Such blocks can be read in and written either from ESPResSo scripts see in the following subsections or from your own C code using the C Interface see section 10 2 1 Writing ESPResSo s global variables Syntax 1 blockfile channel write variable varnamel varname2 2 blockfile channel write variable all Description Variant 1 writes the global variables varnamel1 varname 2 which are known to the setmd command see section 6 1 on page 74 to channel Variant 2 will write all known global variables Note that when the block is read all variables with names listed in the Tcl variable blockfile_variable_blacklist are ignored 10 2 2 Writing Tcl variables Syntax 1 blockfile channel write tclvariable varnamel varname2 2 blockfile channel write tclvariable all 3 blockfile channel write
196. les should most probably not be included the simulation state ESPResSo has no way to distinguish between these variables On the other hand the ESPResSo core has a number of internal variables e g the particle coordinates While most of these are probably good candidates for being included into a checkpoint this is not necessarily so For example when you have particles in your system that have fixed coordinates should these be stored in a checkpoint or not If the system contains mostly fixed particles and only very few moving particles this would increase the memory size of a checkpoint needlessly And what about the interactions in the system or the bonds Should these be stored in a checkpoint or are they generated by the script Another problem with a generic checkpoint would be the control flow of the script In principle the checkpoint would have to store where in the script the checkpointing function was called to be able to return there All this is even further complicated by the fact that ESPResSo is running in parallel Instead in ESPResSo the user has to specify what information needs to be saved to a file to be able to restore the simulation state The blockfile command helps you to do that 10 2 blockfile Using the structured file format ESPResSo uses a standardized ASCII block format to write structured files for anal ysis or storage Basically the file consists of blocks in curled braces which have a single word tit
197. line like the following to the analysis tcl file source file join file dirname info script myanalysis tcl You then need to implement the following essential functions within your new names pace e mbtools analysis myanalysis setup_myanalysis args Typically you would use this function to initialise variables and open files Called by mbtools analysis setup_analysis Arguments are allowed e mbtools analysis myanalysis printav_myanalysis void This function should print results to a file Called by mbtools analysis print_averages Arguments are not allowed e mbtools analysis myanalysis analyze_myanalysis void This function performs the actual analysis and should update the storage and averaging variables Called by mbtools analysis do_analysis Arguments are not allowed e mbtools analysis myanalysis resetav_myanalysis void This function should update averages and reset variables accordingly depending on your requirements Called by mbtool1s analysis reset_averages Arguments are not allowed If any of these functions is not implemented the program will probably crash 14 5 System generation Package for setting up lipid membrane systems in a variety of geometrical shapes 166 14 5 1 Basic commands mbtools system_generation setup_system system_specs ibox1 moltypes e system_specs This is a list of structures called system specifications Each
198. ll 3D mapping tool but this will only rarely be used The full syntax is Syntax observable new needs_radial_profile_specs other_parameters center lt cx gt lt cy gt lt cx gt maxr mazr minz minz maxz mazz rbins rbins phibins phibins zbins zbins 119 Description 9 2 Correlations 9 2 1 Introduction Time correlation functions are ubiquitous in statistical mechanics and molecular simu lations when dynamical properties of many body systems are concerned A prominent example is the velocity autocorrelation function v t v t 7 which is used in the Green Kubo relations In general time correlation functions are of the form C r A t B t 7 9 1 where t is time 7 is the lag time time difference between the measurements of vector observables A and B and is an operator which produces the vector quantity C from A and B The ensemble average is taken over all time origins t Correlation functions describing dynamics of large and complex molecules such as polymers span many orders of magnitude ranging from MD time step up to the total simulation time ESPResSo uses a fast correlation algorithm see section 9 2 6 which enables efficient computation of correlation functions spanning many orders of magnitude in the lag time The generic correlation interface of ESPResSo may process either observables defined in the kernel or data which it reads from an external file or values en
199. llowsphere typeinfo moltypeid hollowsphere sphereparticlelist bondtype natomsfill e sphereparticlelist A list of the particle types for each atom in the hollowsphere The atoms that make up the outer shell must be listed first followed by the atoms that make up the inner filling e bondtype The typeid for bonds linking atoms in the outer shell e natomsfill Number of filler atoms The atom types for these will be obtained from the last natomsfill in the sphereparticlelist 172 Creates a sphere of beads arranged such that they have an approximate spacing of bondl and such that they optimally cover the sphere The optimal covering is obtained using the icover routines which are copyright R H Hardin N J A Sloane and W D Smith 1994 2000 Thus the routine will only work if you have installed icover and if you can successfully run it from the command line in the directory that you started your espresso job These routines are serious overkill so if anybody can think of a nice simple algorithm for generating a covering of the sphere let us know protein typeinfo moltypeid protein particletypelist bondtypelist e particletypelist A list of the particle types for each atom in the protein e bondtypelist A list of bondtypeids Create a protein molecule spanlipid typeinfo moltypeid protein particletypelist bondtypelist e particletypelist A list of the particle types for each atom in the lipid Since
200. lly determine whether an MPI compiler is available If it is it will use it If you specify without mpi or with mpi no then MPI will not be used even if it is available with efence without efence Whether or not to use the electric fence mem ory debugging library Efence is not used by default http freshmeat net projects efence 21 with tcl TCL By default configure will automatically determine which version of Tcl is used If the wrong version is chosen automatically you can specify the name of the library with this option e g tc18 4 with tk TK without tk By default the GUI toolkit Tk is not used by ESPResSo This option can be used to activate Tk and to specify which Tk version to use e g tk8 4 If you only specify with tk and do not give a version number configure will try to automatically deduce the right version with fftw without fftw This can be used to specify whether the FFTW li brary is to be used and which version By default version 3 will be used if it is found otherwise version 2 is used Note that quite a number of central features of ESPResSo require FFTW with cuda path without cuda This switch enables CUDA support path should be the path to the CUDA directory which can be omitted if it is the NVIDIA de fault path i e usr local cuda The variable NVCCFLAGS can be used to define compiler flags for the NVIDIA CUDA compiler nvcc For example NVCCFLAGS
201. lm James L Harden and Gary W Slater Implicit Method for Simulating Electrohydrodynamics of Polyelectrolytes Phys Rev Lett 105 14 SEP 29 2010 doi 10 1103 PhysRevLett 105 148301 R W Hockney and J W Eastwood Computer Simulation Using Particles IOP 1988 F H fling Karl Ulrich Bamberg and Thomas Franosch Anomalous transport resolved in space and time by fluorescence correlation spectroscopy Soft Matter 7 1358 2011 Alan M Horowitz A generalized guided monte carlo algorithm Phys Lett B 268 247 252 1991 219 26 27 28 29 33 34 35 W Humphrey A Dalke and K Schulten VMD Visual molecular dynamics J Mol Graphics 14 33 38 1996 C Junghans and S Poblete A reference implementation of the adaptive resolution scheme in ESPResSo Comp Phys Comm 181 8 1449 1454 2010 C Junghans M Praprotnik and K Kremer Transport properties controlled by a thermostat An extended dissipative particle dynamics thermostat Soft Matter 4 156 2008 Stefan Kesselheim Marcello Sega and Christian Holm Applying to dna translocation Effect of dielectric boundaries Computer Physics Com munications 182 1 33 35 2011 ISSN 0010 4655 doi 10 1016 j cpc 2010 08 014 URL http www sciencedirect com science article pii S001046551000305X ce title Computer Physics Communications Special Edi tion for Conference on Computational Physics Kaohsiung
202. lobal area conservation global variables box_l cell_grid cell_size dpd_gamma dpd_r_cut gamma integ_switch local_box_1 max_cut_bonded 74 max_cut_nonbonded max_cut max_num_cells max_part max_range max_skin min_global_cut 75 min_num_cells n_layers n_nodes n_part_types n_part node_grid npt_p_ext 75 npt_p_inst nptiso_gamma0 nptiso_gammav periodicity piston skin temperature 75 thermo_switch time_step 76 time timings 76 transfer_rate 76 verlet_flag 76 verlet_reuse 76 gyration tensor harmonic bond hat interaction Hertzian interaction hydrodynamic radius of a chain ICCx iccp3m Tcl command icosaeder Tcl command IMD 136 225 imd Tcl command Installation integ_switch global variable integrate Tcl command inter Tcl command 43 Interaction DPD interactions bending force BMHTF bond angle bonded bonded fsi Buckingham Coulomb DAWAANR method Debye Hiickel dihedral Dipolar Directional Lennard Jones DLC method DPD ELC method Electrostatic FENE gaussian Gay Berne Generic Lennard Jones global area conservation harmonic hat hertzian Lennard Jones Lennard Jones cosine local area conservation Maggs method Magnetostatic MDDS method MEMD MMMI1D MMM2D Morse non bonded P3M rigid bond smooth step soft sphere stretching for
203. lobal variable thermo_switch global eer thermostat Tcl command time global variable time unit 9 time_step global variable timings global variable topologies transfer_rate global variable Tunable slip boundary interaction units 9 uwerr Tcl command van Hove autocorrelation function G r t 101 vcf velocities Tcl command verlet_flag global variable verlet_reuse global variable virtual sites volume conservation vsf vtf whitespace 128 writepdb Tcl command writepdbfoldchains Tcl command writepdbfoldtopo Tcl command 228
204. locity of the center of mass of the whole system from every particle s velocity thereby performing a Galilei transform into the reference frame of the center of mass of the system This transformation is useful for example in combination with the DPD thermostat since there a drift in the velocity of the whole system leads to an offset in the reported temperature 88 7 Running the simulation 7 1 integrate Running the simulation Syntax 1 integrate steps 2 integrate set nvt 3 integrate set npt_isotropic Pert piston x y 2 cubic_box Description ESPResSo uses the Velocity Verlet algorithm for the integration of the equations of motion The command integrate with an integer steps as parameter integrates the system for steps time steps Two methods for the integration can be set For an NVT ensemble thermostat and for an NPT isotropic ensemble barostat The current method can be detected with the command integrate set without any parameters The NVT integrator is set without parameters the temperature can be set with the thermostat For the NPT ensemble the parameters that can be added are e Pex The external pressure as float variable This parameter is required e piston The mass of the applied piston as float variable This parameter is required e xyz Three integers to set the box geometry for non cubic boxes This parameter is optional e cubic_box If this optional parameter is added a cubic box is assu
205. ly nptiso_gamma0 double read only nptiso_gammav double read only npt_p_ext double read only Pressure for NPT simulations npt_p_inst double Pressure calculated during an NPT isotropic integration piston double read only Mass off the box when using NPT isotropic integrator periodicity bool 3 Specifies periodicity for the three directions If the feature PARTIAL_PERIODIC is set this variable can be set to 1 1 1 or 0 0 0 at the moment If not it is readonly and gives the default setting 1 1 1 skin double Skin for the Verlet list temperature double read only Temperature of the simulation thermo_switch double read only Internal variable which thermostat to use time double The simulation time 79 time_step double Time step for MD integration timings int Number of samples to time average over transfer_rate int read only Transfer rate for VMD connection You can use this to transfer any integer value to the simulation from VMD verlet_flag bool Indicates whether the Verlet list will be rebuild The program decides this normally automatically based on your actions on the data verlet_reuse double Average number of integration steps the verlet list has been re used 6 2 thermostat Setting up the thermostat The thermostat command is used to change settings of the thermostat The different available thermostats will be described in the following subsections Note that for a si
206. ly used cell system Only the force energy and pressure loops are implemented separately for each cell model as explained above The domain decomposition or link cell algorithm is implemented in ESPResSo such that the cells equal the ESPResSo cells i e each cell is a separate particle list For an example let us assume that the simulation box has size 20 x 20 x 20 and that we assign 2 processors to the simulation Then each processor is responsible for the particles inside a 10 x 20 x 20 box If the maximal interaction range is 1 2 the minimal possible cell size is 1 25 for 8 cells along the first coordinate allowing for a small skin of 0 05 If one chooses only 6 boxes in the first coordinate the skin depth increases to 0 467 In this example we assume that the number of cells in the first coordinate was chosen to be 6 and that the cells are cubic ESPResSo would then organise the cells on each node in a 6 x 12 x 12 cell grid embedded at the centre of a 8 x 14 x 14 grid The additional 182 cells around the cells containing the particles represent the ghost shell in which the information of the ghost particles from the neighbouring nodes is stored Therefore the particle information stored on each node resides in 1568 particle lists of which 864 cells contain particles assigned to the node the rest contain information of particles from other nodes a Classically the link cell algorithm is implemented differently Instead of having sep arat
207. m The various options are equivalent to those described by analyze pressure in 8 1 21 on page 104 It is called a stress tensor but the sign convention follows that of a pressure tensor The stress tensor is calculated by k l k D p a mio yl ee Fij Tij V V where the notation is the same as for analyze pressure in 8 1 21 on page 104 and the superscripts k and correspond to the components in the tensors and vectors Note that the angular velocities of the particles are not included in the calculation of the stress tensor The command is implemented in parallel 8 3 Output format variant 1 pressure total_pressure_tensor ideal ideal_gas_pressure_tensor bond_type pressure_tensor nonbonded_type pressure_tensor coulomb pressure_tensor specifying the pressure tensor the ideal gas pressure tensor the contributions from bonded interactions the contributions from non bonded interactions and the electro static contributions 8 1 23 Local Stress Tensor Syntax analyze local_stress_tensor periodic_x periodic_y periodic_z range_start_x range_start_y range_start_z range_x range_y range_z bins_x bins_y bins_z Description Computes local stress tensors in the system A cuboid is defined starting at the coordi nate range_start_x range_start_y range_start_z and going to the coordinate range_start_a range_a 106 range_start_y range_y range_start_z range_z
208. m Y 9i9g 6 xij Yijs 24 kimiga Pi Nkim apo _ og is given by a 2 y 2 a a y z for x y z 4 0 and g 0 0 0 g 0 0 0 where N id e blPij nkim e7 Prim bala y 2 k l m 40 Mklm The limit exists but differs for three dimensionally periodic systems by some multi ple of the square of the dipole moment from the spherical limit as obtained by the Ewald summation 48 From the physical point of view the Coulomb interaction is replaced by a screened Coulomb interaction with screening length 1 E is then the energy in the limit of infinite screening length But because of the conditional convergence of the elec trostatic sum this is not necessarily the same as the energy of an unscreened system Since the difference to the Ewald methods only depends on the dipole moment of the system the correction can be calculated easily in linear time and can be ignored with respect to accuracy as well as to computation time For one or two dimensionally systems however E E i e the convergence factor approach equals the spherical summation limit of the Ewald sum and MMM1D and MMM2D do not require a dipole correction Starting from this convergence factor approach Strebel constructed a method of com putational order O N log N which is called MMM 51 The favourable scaling is ob tained very much like in the Ewald case by technical tricks in the calculation of the far 210 formula The far
209. m values of the original data are compressed to produce a single data point in the higher compression level Thus the lag time between the neighbouring values in the higher compression level increases by a factor of m while the number of stored values decreases by the same factor and the number of correlation operations at this level reduces by a factor of m 124 Compression 47 Sa za 1 NE p 2 e level oa a 7 oor ae i fa t 1 i 1 aN y v 0 i 0 i 1 i 2 i 3 i 4 i 5 f ee eee i p 2 i p 1 a se nee e api O A a E Tt pt2 gt Wee Sp ES 1 i 0 i 2 i 4 i 6 i p j p 2 s i 2 p 1 Illia SS ao 4 p 1 i ee eee eae Tt 2 p 2 TS ey Y 2p p 2 E 2 i 0 i 4 seer eee i 2p i 2 p 2 i 4 p 1 Figure 9 1 Schematic representation of buffers in the correlator Correlations for lag times 7 2p 4 p 1 are computed at compression level 2 which is created in an analogous manner from level 1 This can continue hierarchically up to an arbitrary level for which enough data is available Due to the hierarchical reduction of the data the algorithm scales as O p log Tmax Thus an additional order of magnitude in Tmax costs just a constant extra effort The speedup is gained at the expense of statistical accuracy The loss of ac
210. med 7 2 change_volume Changing the box volume Syntax 1 change_volume View 2 change_volume Inew x y z xyz Description Changes the volume of either a cubic simulation box to the new volume Vhew or its given x y z xyz extension to the new box length Lnew and isotropically adjusts the particles coordinates as well The function returns the new volume of the deformed simulation box 89 7 3 Stopping particles Use the following functions also see Section 6 9 e kill_particle_motion halts all particles in the current simulation setting their velocities to zero as well as their angular momentum if the feature ROTATION has been compiled in e kill_particle_forces sets all forces on the particles to zero as well as all torques if the feature ROTATION has been compiled in 7 4 velocities Setting the velocities Syntax velocities Umax start pid count N Description Sets the velocities of the particles with particle IDs between pid and pid N to a random vector with a length less than nax and returns the absolute value of the total velocity assigned By default all particles are affected 7 5 invalidate_system Syntax invalidate_system Description Forces a system re init which among others causes the integrator to also update the forces at its beginning instead of re using the values from the previous integration step This is particularly necessary to ensure continuity after settin
211. mentioned in the correct order Command part O bond 33 1 2 3 56 creates a bond related to the angle between the triangles 012 and 123 The particle 0 corresponds to point Al particle 1 to C particle 2 to B and particle 3 to A2 There are two rules that need to be fulfilled e there has to be an edge between particles 1 and 2 e orientation of the triangle 012 that is the normal vector defined as a vector product 01 x 02 must point to the inside of the immersed object Notice that also concave objects can be defined If 09 is larger than 7 then the inner angle is concave 5 4 3 Local area conservation Syntax inter bondid area_force_local Da Kal Description This interaction conserves the area of the triangles in the triangulation This type of interaction is available for closed 3D immersed objects as well as for 2D sheet flowing in the 3D flow The deviation of the triangle surface Sapo is computed from the triangle surface in the resting shape AS4pc SABC Sa spc The area constraint assigns the following shrinking expanding force to every vertex ASABC Fa A kal Sine WA 5 24 where k is the area constraint coefficient and wa is the unit vector pointing from the centroid of triangle ABC to the vertex A Similarly the analogical forces are assigned to B and C This interaction is symmetric therefore after defining the interaction inter 44 area_force_local 0 02 4 0 the following commands ar
212. mulation of the NPT ensemble you need to use a standard thermostat for the particle velocities e g Langevin or DPD and a thermostat for the box geometry e g the isotropic NPT thermostat You may combine different thermostats at your own risk by turning them on one by one Note that there is only one temperature for all thermostats Syntax 1 thermostat 2 thermostat off 3 thermostat parameters Description Variant 1 returns the thermostat parameters A Tcl list is given containing all the parameters needed to set the specific thermostat exactly the same as the input command line without the preceding thermostat Variant 2 turns off all thermostats and sets all thermostat variables to zero Setting temperature to zero also affects the way in which electrostatics are handled see also Section 5 7 Variant 3 sets up one of the thermostats described below 6 2 1 Langevin thermostat Syntax thermostat langevin temperature gamma Description The Langevin thermostat consists of a friction and noise term coupled via the fluctuation dissipation theorem The friction term is a function of the particle velocities For a more detailed explanation refer to If the feature ROTATION is compiled in the rotational degrees of freedom are also coupled to the thermostat 76 Using the Langevin thermostat it is posible to set a temperature and a friction coef ficient for every particle individually Consult the part c
213. n approach requires the VIRTUAL_SITES_RELATIVE feature e The bind at point of collision approach cannot handle collisions between virtual sites 6 8 Catalytic Reactions With the help of the feature CATALYTIC_REACTIONS one can define three particle types to act as reactant catalyzer and product Using these reaction categories we model the following chemical reaction system rt pr 6 3 ct rt pr 6 4 where the first line indicates that there is a reversible chemical reaction that converts the reactant particles rt into product pr particles leading to an equilibrium state The second line indicates that in the presence of a catalyst ct the forward reaction pathway is favored i e conversion of reactants into products The equilibrium reaction is described by the equilibrium constant keq A pr Ka Gt E 6 5 with rt and pr the reactant and product concentration and keq the forward and backward reaction rate constants respectively The rate constants that specify the change in concentration for the equilibrium and catalytic reaction are given by an Keq PT keq rt 6 6 Merl kairt hea Ler 6 7 dirt dlpt el PA ld 6 8 85 respectively In the current ESPResSo implementation we assume keq keq k eq and therefore Keq 1 The user can specify k_eq gt 0 and k_ct k gt 0 The former rate constant is applied to all reactant and product pa
214. n initially some particles are placed on top of each other Using the Tcl interface one can simulate the randomly set up system with capped forces interactively check whether it is safe to remove the cap and switch on the full interactions and then perform the actual productive simulation 1 2 Available simulation methods e Ensembles NVE NVT NpT Algorithms for charged systems P3M for fully periodic systems ELC and MMM family of algorithms for charged systems with non periodic boundary conditions MEMD Maggs algorithm Hydrodynamics DPD as a thermostat Lattice Boltzmann Non equilibrium MD to simulate shear flow e Parallel tempering Metadynamics e Rigid bodies via virtual sites e AdResS 1 3 Basic program structure As already mentioned ESPResSo consists of two components The simulation engine is written in C for the sake of computational efficiency The steering or control level is interfaced to the kernel via an interpreter of the Tcl scripting language The kernel performs all computationally demanding tasks Before all integration of Newton s equations of motion including calculation of energies and forces It also takes care of internal organization of data storing the data about particles communication between different processors or cells of the cell system The kernel is modularized so that basic functions are accessed via a set of well defined lean interfaces h
215. n the C core of ESPResSo during a call to the function integrate while they are set up and their results are collected from the Tcl level These observables are more complex to im plement and offer less flexibility while the are significantly faster and more memory efficient and they can be set up to be computed every few timesteps The observables in this class are described in chapter 9 The class of Tcl level analysis functions is mainly controlled via the analyze command It has two main uses Calculation of observables analyze observable and definition and analysis of topologies in the system analyze topologycommand In addition ESPResSo offers the command uwerr see section 8 4 for computing statistical errors in time series 8 1 Available observables The command analyze provides online calculation of local and global observables 8 1 1 Minimal distances between particles Syntax 1 analyze mindist type_list_a type_list_b 2 analyze distto pid 3 analyze distto z y 2 Description Variant 1 returns the minimal distance between two particles in the system If the type lists are given then the minimal distance between particles of only those types is determined distto returns the minimal distance of all particles to particle pid variant 2 or to the coordinates x y z Variant 3 97 8 1 2 Particles in the neighbourhood Syntax 1 analyze nbhood pid r_catch 2 analyze nbhood z y z ratch Descr
216. natively a type number can be given Exporting the positions of all particles but in separate files might make sense if one wants to distinguish the different particle types in the visualization i e by color or size 10 5 Writing PDB PSF files The PDB Brookhaven Protein DataBase format is a widely used format for describing atomistic configurations PSF is a format that is used to describe the topology of a PDB file When visualizing your system with VMD it is recommended to use the VTF format instead see section 10 3 as it was specifically designed for visualizations with VMD In contrast to the PDB PSF formats in VTF files it is possible to specify the VDW radii of the particles to have a varying simulation box size etc 10 5 1 writepsf Writing the topology Syntax writepsf file molecule Np MPC NcoI N S N S Description Writes the current topology to the file file here file is not a channel since additional information cannot be written anyway Np MPC and so on are parameters describing a system consisting of equally long charged polymers counterions and salt This infor mation is used to set the residue name and can be used to color the atoms in VMD If you specify molecule the residue name is taken from the molecule identity of the particle Of course different kinds of topologies can also be handled by modified versions of writepsf 10 5 2 writepdb Writing the coordinates Syntax 1 writepdb file 2 w
217. nd the point of contact or you can use a virtual bond to prevent additional force contributions at the expense of RATTLE see 5 3 4 The parameters d and bond1 are the same as for the bind_centers method bond2 determines the type of the bond created between the virtual sites if applicable and can be either a pair or a triple angle bond Tf it is a pair bond it connects the two virtual particles otherwise it constraints the angle between the two real particles around the virtual ones type denotes the particle type of the virtual sites created at the point of collision if applicable Be sure not to define a short ranged interaction for this particle type as two particles will be generated in the same place The code can throw an exception background error in case two particles collide for the first time if the exception keyword is added to the invocation In conjunc tion with the catch command of Tcl this can be used to intercept the collision 84 if catch integrate 0 err foreach exception lrange err 2 end if lrange exception 0 2 collision between particles set i lindex exception 3 set j lindex exception 5 puts particles i and j collided The following limitations currently apply for the collision detection e The method is currently limited to simulations with a single cpu e No distinction is currently made between different particle types e The bind at point of collisio
218. ne those parameters that are not set to a predetermined value using the optional parameters of the tuning command The two tuning methods follow different methods for determining the optimal param eters While the tune version tests different values on a grid in the parameter space the tunev2 version uses a bisection to determine the optimal parameters In general for small systems the tune version is faster while for large systems tunev2 is faster The results of tunev2 are always at least as good as the parameters from the tune version and normally the obtained accuracy is much closer to the desired value 62 During execution the tuning routines report the tested parameter sets the correspond ing k space and real space errors and the timings needed for force calculations the setmd variable timings controls the number of test force calculations Since the error depends on Teut box_l and abox_l the output is given in these units Note that the previous setting of reut cao and mesh will be remembered If you want to retune your electrostatics e g after a major system change you should use inter coulomb lg p3m tune accuracy acc r_cut O mesh O cao O Additional P3M parameters Syntax inter coulomb epsilon metallic epsilon n_interpol points mesh_off zoff yoff zoff Description Once P3M algorithm has been set up it is possible to set some additional P3M param eters with this command The different parameters have the
219. nel channel channel channel write configs write start tag write end write tag arg read start read toend read particles interactions bonds variable seed random bitrandom configs blockfile channel read auto writevs channelld short verbose radius radii auto typedesc typedesc ignore_charges writevc channelld short verbose folded absolute pids pids all userdata userdata vtfpid pid writevtk filename all type writepsf file molecule Np MPC NcoI N S NnS writepdb file writepdbfoldchains file chain_start n_chains chain_length box_l writepdbfoldtopo file shift imd connect port imd positions unfolded fold_chains imd listen seconds imd disconnect prepare_vmd_connection filename wait start constraints lbfluid gpu tau lb_timestep ext_force ist friction gamma gamma_odd gamma_odd Required features LB LB_GPU lagrid agrid 1 or or or visc viscosity jo dens density 112 pulk_vise bulk_viscosity l or 2 pa paa gamma_even gamma_even lbfluid print_interpolated_velocity x y z lbfluid save_ascii_checkpoint filename lbfluid save_binary_checkpoint filename lbfluid load_ascii_checkpoint filename lbfluid load_ binary_checkpoint filename thermostat lb Required features or T liB 2LB GPU lbnode z y z print set args Required features LB lbfluid print vtk pr
220. ng from particles p2 to p1 return p1 pos p2 pos bond_vec_min lt p1 gt lt p2 gt box Calculate bond vector pointing from particles p2 to p1 return MinimumImage p1 pos p2 pos bond_length lt p1 gt lt p2 gt Calculate bond length between particles p and p2 bond_length_min lt p1 gt lt p2 gt box Calculate minimum image bond length between particles pi and p2 bond_angle lt p1 gt lt p2 gt lt p3 gt type Calculate bond angle between particles p1 p2 and p3 If type is r the return value is in radiant If it is d the return value is in degree The default for type is r bond_dihedral lt p1 gt lt p2 gt lt p3 gt lt p4 gt typel 2 2 r Calculate bond dihedral between particles p1 p2 p3 and p4 If type is the return value is in radiant If it is d the return value is in degree The default for type is r part_at_dist lt p gt lt dist gt return position of a new particle at distance dist from p with random orien tation part_at_angle lt p1 gt lt p2 gt lt phi gt len return position of a new particle at distance len default 1 0 from p2 which builds a bond angle phi for p1 p2 p new part_at_dihedral lt p1 gt lt p2 gt lt p3 gt lt theta gt phi len return position of a new particle at distance len default 1 0 from p3 which builds a bond angle phi default random for p2 p3 p new and a dihedral angle theta for p1
221. ning the number of charges in the wire Update y field De 1 wire Ey By ar D 9 3 Subtract the charge qwire from the each charge on the sites of Yywire plane Update x field r y 1 vertex This scheme is repeated until the fields are completely relaxed i e the energy is minimized During repetition the spatial dimensions are permutated to avoid a drift in one direction first step third step Figure D 2 Recursive solution of Gauss law D 4 Time integrator For the time discretization we have adopted the elegant solution which was found by Rottler and Maggs and allows to conserve both time reversibility and phase space volume conservation 1 Update the particle momenta by half a time step 206 2 Update the B field by half a time step 3 Update the particle positions in x direction by half a time step 4 Update the electric field in x direction by half a time step 5 Update the particle positions in y direction by half a time step 6 Update the electric field in y direction by half a time step 7 Update the particle positions in z direction by half a time step 8 Update the electric field in z direction by a full time step 9 Update the particle positions in z direction by half a time step E 10 Update the electric field in y direction by half a time step El 11 Update the particle positions in y direction by half a time step E 12 Update the electr
222. nner volume for example blood cells magnetic beads capsules The boundary of an object is covered with triangular mesh The vertices of the mesh are put into ESPResSo as particles The edges of the mesh will define elastic forces keeping the shape of the object The movement of object will be achieved by adding forces to the mesh points Modelled elastic or rigid objects are immersed in the LB fluid flow The fluid interacts with an elastic object resulting in its deformation this imediately generates forces acting back on the fluid The aim is to describe the immersed object using the notion of particles and to create bonds between these particles representing elastic or rigid forces The objects are composed of a membrane encapsulating the fluid inside the object For now the inside fluid must have the same density and viscosity as the outside fluid The object is represented by its membrane boundary that is discretized using a tri angulation Such a triangulation defines interacting particles distributed on the surface of the immersed object 18 e between two particles corresponding to the edges in the triangulation modelling the stretching of the membrane e between three particles corresponding to the triangles of the triangulation local area or local surface preservation of the membrane e between four particles corresponding to two triangles from the triangulation shar ing a common edge bending of the membrane T
223. non virtual particle O denotes the orientation of the vector n with respect to the non virtual particle s body fixed frame and d the distance between virtual and non virtual particle In words The virtual site is placed at a fixed distance from the non virtual particle When the non virtual particle rotates the virtual sites rotates on an orbit around the non virtual particle s center To use this implementation of virtual sites activate the feature VIRTUAL_SITES_ RELATIVE in myconfig h see sec 3 4 To set up a virtual site 1 Place the particle to which the virtual site should be related It needs to be in the center of mass of the rigid arrangement of particles you create Let its particle id be n 2 Place a particle at the desired relative position make it virtual and relate it to the first particle part v pos pos virtual 1 vs_auto_relate n 3 Repeat the previous step with more virtual sites if desired 4 To update the positions of all virtual sites call integrate 0 Please note e The relative position of the virtual site is defined by its distance from the non virtual particle the id of the non virtual particle and a quaternion which defines the vector from non virtual particle to virtual site in the non virtual particle s body fixed frame The first two are saved in the virtual site s vs_relative attribute while the latter is saved in the quaternion attribute Take care not to overwrite these a
224. ntaining the user data for every particle The user data is appended to the coordinate line and can be read into VMD via the VMD plugin VIFTools The default is to provide no userdata Example userdata red blue green 3 3 vtfpid Translating ESPResSo particles ids to VMD particle ids Syntax vtfpid pid Description If pid is the id of a particle as used in ESPResSo this command returns the atom id used in the VIF VSF or VCF formats 10 Thi Par 4 writevtk Particle Visualization in paraview s feature allows to export the particle positions in a paraview compatible VTK file aview is a powerful and easy to use open source visualization program for scientific data Since ESPResSo can export the lattice Boltzmann velocity field in the VTK format as well and paraview allows to visualize particles with glyphs and vector fields with stream lines glyphs contourplots etc one can use it so completely visualize a pled lattice Boltzmann MD simulation It can also create videos without much effort if one exports data of individual time steps into separate files with filenames including cou a running index data_0 vtk data_1 vtk Syntax writevtk filename all type http www paraview org 134 Description Arguments e filename Name of the file to export the particle positions into e all type Specifies which particle type should be exportet The default is all Alter
225. ntains the output from your simulation To view the results cd to the simplebi layer directory and look at the contents The directory contains many files including e The configurations generated during warmup warm gz e pdb files corresponding to warmup configurations warm vmd gz The configurations generated during the main run simplebilayer gz e pdb files corresponding to main run configs simplebilayer vmd gz e The most recently generated checkpoint file checkpoint latest gz e A directory containing the second most recent checkpoint file checkpoint_bak e A file containing the topology of the system simplebilayer top The original parameter file with which you ran the simulation simplebilayer tcl 159 A original parameter file with which you ran the simulation simplebilayer tcl Files containing analysis output for example time_vs_box1_tmp e Force and energy tables tab e VMD script for visualising the warmup warm_animation script e VMD script for visualising the main trajectory vmd_animation script To visualise your results using the vmd scripts you need to make sure that you have vmd installed properly and that you have the special vmd procedures used by the espresso team i e support for the loadseries command Then you can visualise by typing vmd e vmd_animation script 14 3 The main tcl script The main tcl file provided in the examples scripts directory is a relatively complete
226. nteraction is again symmetric After the definition of the interaction by 58 inter 22 volume_force 65 3 3 0 the order of vertices is crucial By the following command the bonds are defined part O bond 22 1 2 Triangle 012 must have correct orientation that is the normal vector defined by a vector product 01 x 02 The orientation must point inside the immersed object 5 5 Bond angle interactions Syntax 1 inter bondid angle_harmonic K pp 2 inter bondid angle_cosine K do 3 inter bondid angle_cossquare K do Required features BOND_ANGLE Description This creates a bond type with identificator bondid with an angle dependent potential This potential is defined between three particles The particle for which the bond is created is the central particle and the angle between the vectors from this particle to the two others determines the interaction K is the bending constant and the optional parameter y is the equilibirum bond angle in radian ranging from 0 to m If this param eter is not given it defaults to py 7 which corresponds to a stretched configuration For example for a bond defined by part p_2 bond 4 p_1 p_3 the minimal energy configurations are the following inter 4 angle_type 1 0 PI inter 4 angle_type 1 0 expr PI 2 P3 o P1 P2 P3 P1 P2 For the potential acting between the three particles three variants are possible e Harmonic bond angle potential 1 A classical harmonic po
227. of VMD which cannot handle multiple atom sections for the same atom However charges are written out per atom so that VMD forgets about atom types and radii If you don t need the charges you can switch them off by this flat making the other settings work again 10 3 2 writevcf Writing the coordinates Syntax writevc channelld short verbose folded absolute pids pids all userdata userdata Description Writes a coordinate or timestep block that contains all coordinates of the system s particles to the channel given by channelld channelld must be an identifier for an open channel such as the return value of an invocation of open Arguments e short verbose Specify whether the output is in a human readable but somewhat longer format verbose or in a more compact form short The default is verbose 133 10 folded absolute Specify whether the particle positions are written in absolute coordinates absolute or folded into the central image of a periodic system folded The default is absolute pids pids all Specify the coordinates of which particles should be writ ten If all is used all coordinates will be written in the ordered timestep format Otherwise pids has to be a Tcl list specifying the pids of the particles The default is all Example pids 0 23 42 luserdata userdata Specify arbitrary user data for the particles userdata has to be a Tcl list co
228. omb or magnetic contributions only Variants 3 and 4 return the energy contributions of the bonded resp non bonded interactions Output format variant 1 energy value kinetic value interaction value 8 1 21 Pressure Syntax 1 analyze pressure 2 analyze pressure total 3 analyze pressure totals ideal coulomb tot_nonbonded_inter tot_nonbonded_intra 4 analyze pressure bonded bondid 5 analyze pressure nonbonded typeidl typeid2 6 analyze pressure nonbonded_intra typeid 7 analyze pressure nonbonded_inter typeid Description Computes the pressure and its contributions in the system Variant 1 returns all the contributions to the total pressure Variant 2 will return the total pressure only Variants 3 4 and 5 return the corresponding contributions to the total pressure The pressure is calculated if there are no electrostatic interactions by 2E kinetic 4 2 gt Pega 8 1 Vf 3V where f 3 is the number of translational degrees of freedom of each particle V is the volume of the system Exinetic is the kinetic energy Fi the force between particles i and 104 j and rj is the distance between them The kinetic energy divided by the degrees of freedom is e e 5 demir 8 2 2 Note that Equation can only be applied to pair potentials and central forces Description of how contributions from other interactions are calculated is beyond the scope of this manual T
229. oment of the specified group of particles y 5 qr Required feature ELECTROSTATICS interacts_with particle_specifications1 particle_specifications2 cutoff For each particle belonging to particle_specifications1 the observable is unity if a neighbour of a type from particle_specifications2 is found within the distance defined by the cutoff If no such neighbour is found the observable is zero The observable has one dimension per each particle of particle_specifications 1 density_profile particle_specifications profile_specifications Compute the density profile within the specified cube For profile specifications see section 9 1 7 lb_velocity_profile particle_specifications profile_specifications Compute the Lattice Boltzmann velocity profile within the specified cube For profile specifications see section flux_density_profile particle_specifications profile_specifications Compute the flux density within the specified cube For profile specifications see section 9 1 7 117 e radial_density_profile Compute the density profile in cylindrical coordinates For profile specifications see section e radial_flux_density_profile Compute the flux density profile in cylindrical coordinates For profile specifica tions see section e l1b_radial_velocity_profile Compute the Lattice Boltzmann velocity profile in cylindrical coordinates For profile specifications see section e tclcommand dimQ command An arbitrary Tcl func
230. ommands reference chapter 4 for information on how to achieve this 6 2 2 GHMC thermostat Syntax thermostat ghmc temperature n_md phi no_flip flip random_flip no_scale scale Description ESPResSo implements Generalized Hybrid Monte Carlo GHMC as a thermostat GHMC is a simulation method for sampling the canonical ensemble 36 The method consists of MC cycles that combine a few constant energy MD steps specified by n_md followed by a Metropolis criterion for their acceptance Prior to integration the particles momenta are mixed with momenta sampled from the appropriate Boltzmann distribution Given the particles momenta p from the last jt GHMC cycle the new momenta are generated by p cos p sin where is a noise vector of random Gaussian variables with zero mean and variance 1 temperature see 25 for more details The momenta mixing parameter cos p corresponds to phi in the implementation In case the MD step is rejected the particles momenta may be flipped This is specified by setting the no_flip flip option for the random_flip option half of the rejected MD steps randomly result in momenta flip The default for momenta flip is no_flip The noise vector s variance van be tuned to exactly 1 temperature by specifying the scale option The default for temperature scaling is no_scale 6 2 3 Dissipative Particle Dynamics DPD ESPResSo implements Dissipative Particle Dynamics
231. on of the thermostat and its random numbers leading to slightly different results compared to the uninterrupted run see The invalidate system command for details The C implementation is t_random 11 2 2 The bit_random command e Without further arguments bit_random returns a random double between 0 and 1 using the R250 generator XOR ing a table of 250 linear independent integers e bit_random seed returns a tcl list with the seeds of the random number generators on each of the n_nodes nodes while bit_random seed lt seed 0 gt lt seed n_nodes 1 gt sets those seeds to the new values respectively re initialising the random num ber generators on each node Note that this is automatically done on invoking Espresso however due to that your simulation will always start with the same ran dom sequence on any node unless you use this tcl command to reset the sequences seeds e Since internally the random number generators random sequences are not based on mere seeds but an array of 250 linear independent integers whose bits are used as matrix elements which are XOR ed to recover the exact state of the random 145 number generators at a given time during the simulation run e g for saving a checkpoint requires knowledge of all these values They can be accessed by bit_random stat which returns a tcl list with all status informations for any node e g 8 nodes gt approx 2016 parameters To overwri
232. onal external force is applied to the particle ext_torque x y z An additional external torque is applied to the particle exclude pid2 Restrictions pid2 must be an existing particle Between the current particle an the exclusion partner s no nonbonded interactions are calculated Note that unlike bonds exclusions are stored with both partners Therefore this command adds the defined exclusions to both partners exclude delete pid2 Searches for the given exclusion and deletes it Again deletes the exclusion with both partners mass mass Sets the mass of this particle to mass If not set all particles have a mass of 1 in reduced units dipm moment Sets the dipol moment of this particle to moment dip dx dy dz Sets the orientation of the dipole axis to dx dy dz 27 e virtual v Declares the particles as virtual 1 or non virtual 0 default Please read chapter 4 4 before using virtual sites e vs_auto_relate_to pid Automatically relates a virtual site to a non virtual particle for the relative implementation of virtual sites pid is the id of the particle to which the virtual site should be related e vs_relative pid distance Allows for manual access to the attributes of virtual sites in the relative implementation pid denotes the id of the particle to which this virtual site is related and distance the distance between non virtual and virtual particle e temp T If used in combination with
233. ons one by one In this construct only two commands are not standard Tcl commands the random number generator t_random and the part command which is used to specify particle properties here the position the charge q and the type In ESPResSo the particle type is just an integer number which allows to group particles it does not imply any physical parameters Here we use it to tag the charges positive charges have type 0 negative charges have type 1 Now we define the ensemble that we will be simulating This is done using the thermostat command We also set some integration scheme parameters setmd time_step 0 01 setmd skin 0 4 set temp 1 set gamma 1 thermostat langevin temp gamma This switches on the Langevin thermostat for the NVT ensemble with temperature temp and friction gamma The skin depth skin is a parameter for the link cell system which tunes its performance but cannot be discussed here Before we can really start the simulation we have to specify the interactions between our particles We use a simple purely repulsive Lennard Jones interaction to model the hard core repulsion 21 and the charges interact via the Coulomb potential set sig 1 0 set cut expr 1 12246 sig set eps 1 0 set shift expr 0 25 eps inter 0 O lennard jones eps sig cut shift O inter 1 0 lennard jones eps sig cut shift O inter 1 1 lennard jones eps sig cut shift O inter coulomb 10 0 p3m tunev2 accuracy le 3 mesh 32 The fi
234. onstraint 4 3 3 Getting the currently defined constraints Syntax constraint num 38 Description Prints out all constraint information If num is specified only this constraint is displayed otherwise all constraints will be printed 4 4 Virtual sites Virtual sites are particles the positions and velocities of which are not obtained by integrating an equation of motion Rather their coordinates are obtained from the position and orientation of one or more other particles In this way rigid arrangements of particles can be constructed and a particle can be placed in the center of mass of a set of other particles Virtual sites can interact with other particles in the system by means of interactions Forces are added to them according to their respective particle type Before the next integration step the forces accumulated on a virtual site are distributed back to those particles from which the virtual site was derived There are two distinct types of virtual sites decribed in the following 4 4 1 Virtual sites in the center of mass of a molecule To activate this implementation enable the feature VIRTUAL_SITES_COM in myconfig h sec 3 4 Virtual sites are then placed in the center of mass of a set of particles as defined below Their velocity will also be that of the center of mass Forces accumulating on the virtual sites are distributed back to the particles which form the molecule To place a virtual site at th
235. operty filename lbboundary shape shape_args velocity vz vy vz Required features LB_BOUNDARIES pun N Ko ma Ko 2 m N Ko m w j m w w N wo w w A ER w a E w al E w al O 3 pue w a ew Kio ew Ko pue j aS ma A 00 19 al lbboundary force boundary Required features LB_BOUNDARIES lbfluid cpu lbfluid gpu Required features 113 LB_GPU setmd mu_E pk uE pE Required features LB LB_ELECTROHYDRODYNAMICS countBonds particle ist findPropPos particle ropertyist property findBondPos particle roperty ist timeStamp path prefix postfix suffix gt Ko m Ko 4 E al o g g Afafa fa o fl fo 196 B Features This chapter describes the features that can be activated in ESPResSo Even if possible it is not recommended to activate all features because this will negatively effect ESPResSo s performance Features can be activated in the configuration header myconfig h see section 3 4 on page 24 Too activate FEATURE add the following line to the header file define FEATURE B 1 General features e PARTIAL_PERIODIC By default all coordinates in ESPResSo are periodic With PARTIAL_PERIODIC turned on the ESPResSo global variable periodic see sec tion 6 1 on page 74 controls the periodicity of the individual coordinates Note that this slows the integrator down by
236. or Required features ELECTROSTATICS Description MMMID coulomb method for systems with periodicity 0 0 1 Needs the nsquared cell system see section 6 4 on page 81 The first form sets parameters manually The switch radius determines at which xy distance the force calculation switches from the near to the far formula If the Bessel cutoff is not explicitly given it is determined from the maximal pairwise error otherwise this error only counts for the near formula The second tuning form just takes the maximal pairwise error and tries out a lot of switching 64 radii to find out the fastest one If this takes too long you can change the value of the setmd variable timings which controls the number of test force calculations For details on the MMM family of algorithms refer to appendix E on page 210 5 7 5 Maxwell Equation Molecular Dynamics MEMD Syntax inter coulomb lg memd f_mass mesh epsilon Required features ELECTROSTATICS Description This is an implementation of the instantaneous 1 r Coulomb interaction A 5 35 r as the potential of mean force between charges which are dynamically coupled to a local electromagnetic field The algorithm currently works with the following constraints e cellsystem has to be domain decomposition but without Verlet lists e system has to be periodic in three dimensions Arguments e f_mass is the mass of the field degree of freedom and equals to the square root of th
237. p a directory for simulation output It copies forcetables and the parameter file to the directory after creating it if necessary mbtools utils read_startfile file e file Complete path of the file to be read Should be an espresso blockfile Read in particle configuration from an existing file or simulation snapshot mbtools utils read_checkpoint dir e dir Directory containing the checkpoint file which must be called checkpoint latest gz Read in a checkpoint and check for success Warn if the checkpoint does not exist mbtools utils read_topology file e file Complete path of the file that contains the topology information Read in the topology from a file and then execute the analyze set topo_part_sync command of ESPResSo mbtools utils set_topology topo e topo A valid topology Set the given topology and then execute the analyze set topo_part_sync command of ESPResSo mbtools utils set_bonded_interactions bonded_parms e bonded arms A list of bonded interactions Each element of this list should contain all the appropriate arguments in their correct order for a particular call to the espresso inter command See the espresso inter command for a list of possible bonded interactions and correct syntax 174 Set all the bonded interactions mbtools utils set_nb_interactions nb_parms e nb_parms A list of interactions Each element of this list should contain all the appropriate ar
238. p2 p3 p new e INTERACTION RELATED Help functions related to interactions implemented in ESPResSo 142 calc_lj_shift lt lj_sigma gt lt lj_cutoff gt returns the value needed to shift the Lennard Jones potential to zero at the cutoff e VECTOR OPERATIONS A vector v is a tcl list of numbers with an arbitrary length Some functions are provided only for three dimensional vectors corresponding functions contain 3d at the end of the name veclen lt v gt return the length of a vector veclensqr lt v gt return the length of a vector squared vecadd lt a gt lt b gt add vector a to vector b return a b vecsub lt a gt lt b gt subtract vector b from vector a return a b vecscale lt s gt lt v gt scale vector v with factor s return s v vecdot_product lt a gt lt b gt calculate dot product of vectors a and b return a b veccross_product3d lt a gt lt b gt calculate the cross product of vectors a and b return a x b vecnorm lt v gt len normalize a vector to length len default 1 0 unitvec lt p1 gt lt p2 gt return unit vector pointing from position p1 to position p2 orthovec3d lt v gt len return orthogonal vector to v with length len default 1 0 This vector does not have a random orientation in the plane perpendicular to v create_dihedral_vec lt v1 gt lt v2 gt lt theta gt phi len create last vector of a dihedral v1 v2 re
239. pon this call ESPResSo allocates the necessary amount of memory and returns an integer id which will be used later to refer to the observable The parameter name and further arguments have to correspond to one of the observables described below Available observables Currently the following observables are implemented Particle specifications see sec tion below define a group of particles from which the observable should be calcu lated They are generic to all observables and are described after the list of observables particle_positions particle_specifications Positions of the particles in the format 11 y1 21 T2 Y2 22 Un Yn Zn The particles are ordered ascending according to their ids particle_velocities particle_specifications Velocities of the particles in the format ut UL UL US vs vs ut up uz The particles are ordered ascending according to their ids particle_forces particle_specifications Forces on the particles in the format fi fe FE f3 ff 8 JZ fi f The particles are ordered ascending according to their ids particle_angular_momentum particle_specifications Angular momenta omega of the particles in the format wt wi wi wi wh wi we wi we The particles are ordered ascending according to their ids com_position particle_specifications blocked size Position of the centre of mass If blocked size is specified the particles are subdi vided into blocks of size si
240. r and in the center of the system the error distribution is highly non homogenous This is unproblematic as long as the maximal error is really much smaller than the thermal energy However one cannot interpret the error simply as an additional error source le 0 01 0 001 z oe ar j Eze cts gu E E pee A AR 0 0001 le 05 0 Figure E 1 Error distribution of the ELC method Figure shows the error distribution of the ELC method for a gap size of 10 of the total system height For MMM2D and MMMID the error distribution is less ho mogenous however also here it is always better to have some extra precision especially since it is computationally cheap 217 F Bibliography Andersen Molecular Dynamics Simulations At Constant Pressure And Or Tem perature J Chem Phys 72 4 2384 2393 1980 ISSN 0021 9606 Hans C Andersen Rattle A velocity version of the shake algorithm for molecular dynamics calculations J Comp Phys 51 24 34 1983 4 2 1 A Arnold and C Holm MMMID A method for calculating electrostatic inter actions in 1D periodic geometries J Chem Phys 123 12 144103 2005 5 7 4 Axel Arnold and Christian Holm A novel method for calculating electrostatic interactions in 2D periodic slab geometries Chem Phys Lett 354 324 330 2002 Axel Arnold and Christian Holm MMM2D A fast and accurate summation methodlimnb for electrostatic interactions in
241. r simulation and this user guide describes how to use this tool However it should be borne in mind that being able to operate a tool is not sufficient to obtain physically meaningful results It is always the responsibility of the user to understand the principles behind the model simulation and analysis methods he is using ESPResSo will not do that for you It is expected that the users of ESPResSo and readers of this user guide have a thorough understanding of simulation methods and algorithms they are planning to use They should have passed a basic course on molecular simulations or read one of the renown textbooks e g 20 It is not necessary to understand everything that is contained in ESPResSo but it is inevitable to understand all methods that you want to use Using the program as a black box without proper understanding of the background will most probably result in wasted user and computer time with no useful output To enable future extensions the functionality of the program is kept as general as possible It is modularized so that extensions to some parts of the program e g im plementing a new potential can be done by modifying or adding only few files leaving most of the code untouched To facilitate the understanding and the extensibility much emphasis is put on read ability of the code Hard coded assembler loops are generally avoided in hope that the overhead in computer time will be more than compensated for by saving m
242. re re2 error_of_re2 Radius of gyration Syntax analyze rg lt rg gt chain_start n_chains chain_length Description Returns the radius of gyration averaged over all chains It is a radius of a sphere which would have the same moment of inertia as the molecule defined as N 1 ae al Ry N y ri E Fem gt 8 4 i 1 where 7 are position vectors of individual particles constituting a molecule and Fm is the position vector of its centre of mass The sum runs over all N particles comprising the molecule For more information see any polymer science book e g 44 If lt rg gt is used the radius of gyration is averaged over all stored configurations see section 8 3 on page 112 Output format rg error_of_rg rg2 error_of_rg2 Hydrodynamic radius Syntax analyze rh lt rh gt chain_start n_chains chain_length 108 Description Returns the hydrodynamic radius averaged over all chains If lt rh gt is used the hydro dynamic radius is averaged over all stored configurations see section 8 3 on page 112 The following formula is used for the computation ALE 85 Ra NEGRA l The above mentioned formula is only valid under certain assumptions For more infor mation see Chapter 4 and equation 4 102 in I7 Output format rh error_of_rh Internal distances Syntax analyze internal_dist lt internal_dist gt chain_start n_chains chain_length Description Returns t
243. riage return will be used after the message The mmsg equivalent of puts Designed for printing of simple status or progress mes sages mmsg err namespace string newline e namespace A namespace Typically this should be the current namespace which one can get via namespace current e string The message you want printed e newline yes Set this to anything other than yes and no carriage return will be used after the message Prints error messages and causes program to exit mmsg warn namespace string newline e namespace A namespace Typically this should be the current namespace which one can get via namespace current e string The message you want printed e newline yes Set this to anything other than yes and no carriage return will be used after the message Prints warning messages mmsg debug namespace string newline e namespace A namespace Typically this should be the current namespace which one can get via namespace current e string The message you want printed e newline yes Set this to anything other than yes and no carriage return will be used after the message Prints debug messages 180 14 7 2 Control commands mmsg does several checks before it decides to print a message For any given message type it checks if that message type is allowed It also checks to see if the namespace given as an argument is in the allowable namespaces list Th
244. ritepdbfoldchains file chain_start n chains chain length box_l 3 writepdbfoldtopo file shift Description Variant 1 writes the corresponding particle data Variant 2 writes folded particle data where the folding is performed on chain centers of mass rather than single particles In order to fold in this way the chain topology 135 and box length must be specified Note that this method is outdated Use variant 3 instead Variant 3 writes folded particle data where the folding is performed on chain centers of mass rather than single particles This method uses the internal box length and topology information from espresso If you wish to shift particles prior to folding then supply the optional shift information shift should be a three member tcl list consisting of x y and z shifts respectively and each number should be a floating point ie with decimal point 10 6 Online visualisation with VMD IMD Interactive Molecular Dynamics is the protocol that VMD uses to communicate with a simulation Tcl md implements this protocol to allow online visual analysis of running simulations In IMD the simulation acts as a data server That means that a simulation can provide the possibility of connecting VMD but VMD need not be connected all the time You can watch the simulation just from time to time In the following the setup and usage of IMD is described 10 6 1 imd Using IMD in the script Syntax 1 imd connect port
245. rm clusters space_dist is the distance between two monomers up to which they are consid ered to belong to the same clusters The three parameters may be connected by scaling 102 arguments Make sure that your results are only weakly dependent on the exact choice of your parameters For the algorithm the coordinates stored in partCfg are used The chain itself is defined by the identity first of its first monomer and the chain length length Attention This function is very specific to the problem and might not give useful results for other cases with similar structures 8 1 19 Finding holes Syntax analyze holes typeidprobe mesh_size Description Function for the calculation of the unoccupied volume often also called free volume in a system Details can be found in Schmitz and Muller Plathe 46 It identifies free space in the simulation box via a mesh based cluster algorithm Free space is defined via a probe particle and its interactions with other particles which have to be defined through LJ interactions with the other existing particle types via the inter command before calling this routine A point of the mesh is counted as free space if the distance of the point is larger than LJ_cut LJ_offset to any particle as defined by the LJ interaction parameters between the probe particle type and other particle types How to use this function Define interactions between all or the ones you are interested in particle types in your
246. rst element of the list must be a string with the name of the zone type and subsequent elements will be further information about the zone Available zones are sphere center radius cuboid center L W H Determines whether the point at pos is outside the zone Parameter center should be a tcl list Returns 1 if it is and 0 if it is not mbtools utils calc_com mol e mol The molecule Calculate the center of mass of a molecule mbtools utils centersofmass_bymoltype moltypes e moltypes A list of molecule type ids Determine the center of mass of every molecule whose type matches an item in the list moltypes Returns a nested list where each element in the list is itself a list of centers of mass for a given moltype 14 7 mmsg mmsg is designed to provide a more controlled way of printing messages than the simple puts commands of Tcl It has an ability to turn on or off messages from particular namespaces 179 14 7 1 Basic commands The following commands represent the standard interface for the mmsg package For consistency one should use these instead of a bare puts to standard out mbtools makes extensive use of these commands mmsg send namespace string newline e namespace A namespace Typically this should be the current namespace which one can get via namespace current e string The message you want printed e newline yes Set this to anything other than yes and no car
247. rst three inter commands instruct ESPResSo to use the same purely repulsive Lennard Jones potential for the interaction between all combinations of the two parti cle types 0 and 1 by using different parameters for different combinations one could simulate differently sized particles The last line sets the Bjerrum length to the value 10 and then instructs ESPResSo to use PM for the Coulombic interaction and to try 14 to find suitable parameters for an rms force error below 107 with a fixed mesh size of 32 The mesh is fixed here to speed up the tuning for a real simulation one will also tune this parameter If we want to calculate the temperature of our system from the kinetic energy we need to know the number of the degrees of freedom of the particles In ESPResSo these are usually 3 translational plus 3 rotational degrees of freedom if the feature ROTATION is activated You can get this number in the following way P if regexp ROTATION code_info set deg_free 6 else set deg_free 3 Now we can integrate the system set integ_ steps 200 for set i 0 i lt 20 incr i set temp expr analyze energy kinetic deg_free 2 0 n_part puts t setmd time E analyze energy total T temp integrate integ_steps This code block is the primary simulation loop and runs 20x integ_steps MD steps Every integ_steps time steps the potential electrostatic and kinetic energies are printed out the
248. rthermore you should not use other special features like part fix or constraints inside the top and middle slabs For further reference on how NEMD is implemented into ESPResSo see 49 Variant 1 chooses the momentum exchange method In this method in each step the n_exchange largest positive x components of the velocity in the middle slab are selected and exchanged with the n_erchange largest negative x components of the velocity in the top slab Variant 2 chooses the shear rate method In this method the targetted x component of the mean velocity in the top and middle slabs are given by L target_velocity shearrate a 6 1 80 where L is the simulation box size in z direction During the integration the x component of the mean velocities of the top and middle slabs are measured Then the difference between the mean x velocities and the target x velocities are added to the x component of the velocities of the particles in the respective slabs Variant 3 will turn off NEMD variant 4 will print usage information of the param eters of NEMD Variant 5 will return the velocity profile of the system in x direction mean velocity per slab Variant 6 will return the viscosity of the system that is computed via _ F LaLy where F is the mean force momentum transfer per unit time acting on the slab Lz Ly is the area of the slab and y is the shearrate n 6 2 6 4 cellsystem Setting up the cell system Thi
249. rticle than shield a new attempt of setting up the whole chain is done up to trymax times The default for the mode is RW the default for the shield is 1 0 and the default for trymax is 30000 which is usually enough for PSAW Depending on the length of the chain for the SAW mode trymax has to be increased by several orders of magnitude e charge valency Sets the valency of the charged monomers If the valency of the charged polymers valency is smaller than 10 the charge is assumed to be zero and the types are set to typeidchargea typetdneutral If charge is not set it defaults to 0 0 e distance dehargea Sets the stride between the indices of two charged monomers This defaults defaults to 1 meaning that all monomers in the chain are charged e types typeidneutral typeidchargea Sets the type ids of the neutral and charged monomer types to be used with the part command If only typeidneutra is defined typeidchargea defaults to 1 If the option is omitted both monomer types default to 0 e bond bondid Sets the type number of the bonded interaction to be set up between the monomers This defaults to 0 Any bonded interaction no matter how many bonding partners needed is stored with the second particle in this bond See chapter e angle 0 x y z Allows for setting up helices or planar polymers and theta are the angles between adjacent bonds xz y and z set the position of the second monomer of the first
250. rticles in the system whereas the latter is applied only to the reactant particles in the vicinity of a catalyst particle Reactant particles that have a distance of r or less to at least one catalyzer particle are therefore converted into product particles with rate constant k_eq k_ct The conversion of particles is done stochastically on the basis of the relevant rate constant k gt 0 Per 1 exp kAt 6 9 with Poy the probability of the conversion and At the integration time step If the equilibrium rate constant is not specified it is assumed that k_eq 0 Syntax 0 reaction reactant_type rt catalyzer_type ct product_type pt range r ct_rate k_ct eq_rate k_eq react_once on off 1 reaction off 2 reaction print Required features CATALYTIC_REACTIONS The current implementation also requires the use of verlet lists and domain decomposi tion Description e Variant 0 defines a reaction with particles of type number rt as reactant type ct as catalyzer and type pt as product The catalytic reaction rate constant is given by k_cf and to override the default rate constant for the equilibrium reaction k_eq 0 one can specify it by k_eq By default each reactant particle is checked against each catalyst particle react_once off However when creating smooth surfaces using many catalyst particles it can be desirable to let the reaction rate be independent of the surface density of these particles That is
251. s x y or z e centerofrotation rotate around given point e beadtypelist only analyze beads of given types e thetabins histogram bins in angle theta 8 1 5 Modes Syntax analyze modes2d Description Analyzes the modes of a configuration Requires that a grid is set and that the system contains more than two particles Output are four numbers in the order htrE htiwm Ore Oru 8 1 6 Lipid orientation Syntax 1 analyze get_lipid_orients 2 analyze lipid_orient_order Description 8 1 7 Bilayers Syntax 1 analyze bilayer_set 2 analyze bilayer_density_profile Description 8 1 8 GPB Syntax analyze cell_gpb Manningparameter outercellradius innercellradius accuracy numberofinteractions Description 8 1 9 Get folded positions Syntax analyze get_folded_positions molecule shift zx y z 99 Description Outputs the folded positions of particles Without any parameters the positions of all particles are given folded to the box length The optional parameter molecule ensures that molecules particle groups are kept intact The optional shift parameters can be used to shift the not separated molecules if needed 8 1 10 Vkappa Syntax analyze Vkappa reset read set V 1 Vko avk y Description 8 1 11 Radial distribution function Syntax analyze rdf lt rdf gt part_type_list_a part_type_list_b rmin rmazx rbins Description Returns its parameters and the ra
252. s with dihedral angle theta and bond angle v2 res phi and length len default 1 0 If phi is ommited or set to rnd then phi is assigned a random value between 0 and 2 Pi e TCL LIST OPERATIONS average lt list gt 143 Returns the avarage of the provided list list_add_value lt list gt lt val gt Add val to each element of list flatten lt list gt flattens a nested list list_contains lt list gt lt val gt Checks wether list contains val returns the number of occurences of val in list e REGRESSION LinRegression lt 1 gt where lis a listof pairs of points 1 x1 y1 x2 y2 LinRegression returns the least square linear fit ax b and the standard errors og and op LinRegressionWithSigma lt 1 gt where lis a list of lists of points in the form x1 y1 s1 x2 y2 s2 where s is the standard deviation of y LinRegressionWithSigma returns the least square linear fit ax b the standard errors Ca and op covariance cov a b and x 11 2 1 t_random e Without further arguments t_random returns a random double between 0 and 1 using the ranl random number gener ator from Numerical Recipes e t_random int lt n gt returns a random integer between 0 and n 1 e t_random seed returns a tcl list with the seeds of the random number generators on each of the n_nodes nodes while t_random seed lt seed 0 gt lt seed n_nodes 1 gt sets those seeds to th
253. s are temperatures and call them accordingly It is possible to use multiple processors via TCP IP networking but the number of processors can be smaller than the number of temperatures Arguments e swap specifies the name of the routine calculating the swap probability for a sys tem The routine has to accept three parameters the id of the system to evaluate and two temperatures T and T2 The routine should return a list containing the energy of the system at temperatures T and T s respectively e perform specifies the name of the routine performing the simulation between two swap tries The routine has to accept two parameters the id of the system to propagate and the temperature T at which to run it Return values are ignored e init specifies the name of a routine initializing a system This routine can for example create the particles perform some intial equilibration or open output files The routine has to accept two parameters the id of the system to initialize and its initial temperature T Return values are ignored e R specifies the number of swap trial rounds in each round neighboring temper atures are tried for swapping alternatingly i e with four temperatures The first swap trial round tries to swap 1 2 and 3 4 the second round 2 3 and so on e master the name of the host on which the parallel tempering master node is run ning e port the TCP IP port on which the parallel_tempering master should listen
254. s section deals with the flexible particle data organization of ESPResSo Due to different needs of different algorithms ESPResSo is able to change the organization of the particles in the computer memory according to the needs of the used algorithms For details on the internal organization refer to section 15 1 on page 182 6 4 1 Domain decomposition Syntax cellsystem domain_decomposition no_verlet_list Description This selects the domain decomposition cell scheme using Verlet lists for the calculation of the interactions If you specify no_verlet_list only the domain decomposition is used but not the Verlet lists The domain decomposition cellsystem is the default system and suits most applica tions with short ranged interactions The particles are divided up spatially into small compartments the cells such that the cell size is larger than the maximal interaction range In this case interactions only occur between particles in adjacent cells Since the interaction range should be much smaller than the total system size leaving out all interactions between non adjacent cells can mean a tremendous speed up Moreover since for constant interaction range the number of particles in a cell depends only on the density The number of interactions is therefore of the order N instead of order N if one has to calculate all pair interactions 6 4 2 N squared Syntax cellsystem nsquare 81 Description This selects th
255. s that prints a warning if particles come too close so that the simulation becomes unphysical e OLD_DIHEDRAL Switch the interface of the dihedral potential to its old less flexible form Use this for older scripts that are not yet adapted to the new interface of the dihedral potential If you want to use bond angle potentials see section 5 5 on page 59 you need the followig features e BOND_ANGLE e BOND_ANGLEDIST e BOND_ENDANGLEDIST 200 B 3 Debug messages Finally there are a number of flags for debugging The most important one are ADDITIONAL_CHECKS Enables numerous additional checks which can detect incon sistencies especially in the cell systems This checks are however too slow to be enabled in production runs MEM_DEBUG Enables an internal memory allocation checking system This produces output for each allocation and freeing of a memory chunk and therefore allows to track down memory leaks This works by internally replacing malloc realloc and free The following flags control the debug output of various sections of Espresso You will however understand the output very often only by looking directly at the code COMM_DEBUG Output from the asynchronous communication code EVENT_DEBUG Notifications for event calls i e the on_ functions in initialize c Useful if some module does not correctly respond to changes of e g global vari ables INTEG_DEBUG Integrator output CELL_DEBUG Cellsystem output
256. system and a fictitious particle type Practically one uses the van der Waals radius of the particles plus the size of the probe you want to use as the Lennard Jones cutoff The mesh spacing is the box length divided by the mesh ize Output format n_holes mean_hole_size max_hole_size free_volume_fraction sizes surfaces element_lists A hole is defined as a continuous cluster of mesh elements that belong to the unoccu pied volume Since the function is quite rudimentary it gives back the whole information suitable for further processing on the script level sizes and surfaces are given in number of mesh points which means you have to calculate the actual size via the corresponding volume or surface elements yourself The complete information is given in the element_ lists for each hole The element numbers give the position of a mesh point in the linear representation of the 3D grid coordinates are in the order x y z Attention the algo rithm assumes a cubic box Surface results have not been tested Requires the feature LENNARD_JONES 103 8 1 20 Energies Syntax 1 analyze energy 2 analyze energy total kinetic coulomb magnetic 3 analyze energy bonded bondid 4 analyze energy nonbonded typeidl typeid2 Description Returns the energies of the system Variant 1 returns all the contributions to the total energy Variant 2 returns the numerical value of the total energy or its kinetic or Coul
257. t does not play a role as it only specifies the parameters only when applying the bond using the bond particle the number of involved particles plays a role The number of involved particles and their order if important is nevertheless specified here for completeness 5 3 1 FENE bond Syntax inter bondid fene K Armax ro Description This creates a bond type with identificator bondid with a FENE finite extension nonlin ear expander interaction This is a rubber band like symmetric interaction betweeen two particles with prefactor K maximal stretching Armax and equilibrium bond length ro The bond potential diverges at a particle distance r rg Afmax and r ro Armax It is given by 5 19 Tmax 2 V r 5K Armas In h C 2 52 5 3 2 Harmonic bond Syntax inter bondid harmonic K R reut Description This creates a bond type with identificator bondid with a classical harmonic potential It is a symmetric interaction between two particles The potential is minimal at particle distance r R and the prefactor is K It is given by V r iK r R 5 20 The third optional parameter rey defines a cutoff radius Whenever a harmonic bond gets longer than Teut the bond will be reported as broken and a background error will be raised 5 3 3 Subtracted Lennard Jones bond Syntax inter bondid subt_1j reserved R Description This creates a bond type with identificator bond
258. tat can be used to input output metadynamics information between walkers at regular intervals Warning the information extracted from print_stat con tains the entire history of the simulation while only the last increment of sampling should be communicated between walkers in order to avoid counting the same samples multiple times Details on implementation As of now only two reaction coordinates have been implemented distance and relative_ z Many different reaction coordinates can be set up and it is rather easy to implement new ones See the code in metadynamics h c for further details The bias functions that are applied to the potential of mean force and the biased force are not gaussian function as in many metadynamics codes but so called Lucy functions See for more details These avoid the calculation of exponentials 96 8 Analysis in Tcl ESPResSo has two fundamentally different classes of observables for analyzing the sys tems On the one hand some observables are computed from the Tcl level In that case the observable is measured in the moment that the corresponding Tcl function is called and the results are returned to the Tcl script In general observables in this class should only be computed after a large number of timesteps as switching forth and back between the C and the Tcl level is costly This chapter describes all observables in this class On the other hand some observables are computed and stored i
259. tclvariable reallyall Description These commands will write Tcl global variables to channel Global variables are those declared in the top scope of the Tcl script or those that were explicitly declared global When reading the block all variables with names listed in the Tcl variable blockfile_tclvariable_blacklist are ignored Variant 1 writes the Tcl global variables varnamel varname2 to channel Variant 2 will write all Tcl variables to the file with the exception of the inter nally predefined globals from Tcl tcl_version argv argv0 argc tcl_interactive auto_oldpath errorCode auto_path errorInfo auto_index env tcl_pkgPath 128 tcl_patchLevel tcl_libPath tcl_library and tcl_platform Variant 3 will even write those 10 2 3 Writing particles bonds and interactions Syntax 1 blockfile channel write particles what range all 2 blockfile channel write bonds range 3 blockfile channel write interactions Description Variant 1 writes particle information in a standardized format to channel what can be any list of parameters that can be specified in part part d print except for bonds Note that id and pos will automatically be added if missing range is a Tcl list of ranges which particles to write The keyword all denotes all known particles Variant 2 writes the bond information in a standardized format to channel The involved particles and bond types must exist and be valid Variant 3
260. te here the distance between pid and pidg an iterative biased force pulls the system away from the values of the reaction coordinate most sampled Ultimately the system is driven in such a way that it self diffuses along the reaction coordinate between the two boundaries here dmin and dmax The potential of mean force or free energy profile can be extracted by reading the profile Arguments e pid ID of the first particle involved in the metadynamics scheme e pido ID of the second particle involved in the metadynamics scheme e dmin 2min Minimum value of the reaction coordinate While dmin must be positive it s a distance Zmin can be negative since it s the relative height of pid with respect to pide e dmax Zmax Maximum value of the reaction coordinate e Dheigsn height of the bias function e bwidath width of the bias function e found Strength of the ramping force at the boundaries of the reaction coordinate interval e dhins Zbins number of bins of the reaction coordinate e profile_list Tcl list of a previous metadynamics profile 95 e force_list Tcl list of a previous metadynamics force Details on usage Variant 1 returns the status of the metadynamics routine Variant 2 turns metady namics off default value Variant 3 sets a metadynamics scheme with the reaction coordinate distance which corresponds to the distance between any two particles of the system e g calculate the potential
261. te particles id pos type close f The created files config_ are human readable and look like tclvariable box_1 10 density 0 7 variable box_l 10 0 10 0 10 0 particles id pos type 0 3 51770181433 4 3208975936 5 30529948918 0 1 3 93145531704 6 58506447035 6 95045147034 1 As you can see such a blockfile consists of several Tcl lists which are called blocks and can store any data available from the simulation Reading a configuration is done by the following simple script set f open filename r while blockfile f read auto eof close f The blockfile read auto commands will set the Tcl variables box_1 and density to the values specified in the file when encountering the tclvariable block and set the box dimensions for the simulation when encountering the variable block The particle positions and types of all 216 particles are restored when the particles block is read Note that it is important to have the box dimensions set before reading the particles to avoid problems with the periodic boundary conditions With these configurations we can now investigate the system As an example we will create a second script which calculates the averaged radial distribution functions g r and g _ r The radial distribution function for a the current configuration can be obtained using the analyze command set rdf analyze rdf 0 1 0 9 expr box_1 2 100 set rlist set rdflist foreac
262. te those internally in Espresso e g upon restoring a checkpoint submit the whole list back using bit_random stat lt status list gt with status list being the tcl list mentioned above without any braces Be careful A complete recovery of the current state of the simulation is only possible if you make sure to include a call to The invalidate_system command after you saved the checkpoint tcl checkpoint_set will do this automatically for you because the integration algorithm re uses the old forces calculated in the previous time step if something has changed in the system or if it has just been read from a file the forces are re derived including application of the thermostat and its random numbers leading to slightly different results compared to the uninterrupted run see The invalidate system command for details Note further that the bit wise display of integers as it is used by this random number generator is platform dependent As long as you stay on the same archi tecture this doesn t matter at all however it wouldn t be wise to use a checkpoint including the state of the R250 to restart the simulation on a different platform most likely the integers will have a different bit muster leading to a completely different random matrix So if you re using this random number generator always remain on the same platform 11 3 Checking for features of ESPResSo In an ESPResSo Tcl script you can get information whet
263. tential V do 5 28 Unlike the two following variants this potential has a kink at dy 7 and accordingly a discontinuity in the force and should therefore be used with caution 59 e Cosine bond angle potential 2 V a K 1 cos 0 5 29 Around do this potenial is close to a harmonic one both are 1 2 do in leading order but it is periodic and smooth for all angles e Cosine square bond angle potential 3 V a cos cos do 5 30 This form is used for example in the GROMOS96 force field The potential is 1 8 0 around Ho and therefore much flatter than the two potentials before 5 6 Dihedral interactions Syntax inter bondid dihedral n K p Description This creates a bond type with identificator bondid with a dihedral potential i e a four body potential In the following let the particle for which the bond is created be particle pa and the other bond partners p1 p3 pa in this order i e part po bond bondid p p3 pa Then the dihedral potential is given by V K 1 cos n p 5 31 where n is the multiplicity of the potential number of minimas and can take any integer value typically from 1 to 6 p is a phase parameter and K is the bending constant of the potential is the dihedral angle between the particles defined by the particle quadrupel p p2 p3 and pa i e the angle between the planes defined b
264. ter ks kb kal kag kv Their mathematical formulations have been taken from The mass of the particles and the friction coefficient can be calibrated using the drag coefficients of the ellipsoidal objects These drag coefficients have known analytical values and the mass and friction can be calibrated to fit this values More details about the calibration is in 11 The elastic parameters are specific to the immersed objects They correspond to their physical values More details about their mechanical and biological meaning is presented in specifically for red blood cells However the proper calibration to fit the experimental data has been performed in 11 13 3 Geometry The membrane of the immersed object is triangulated In samples object in fluid you can find an example using deformable objects in the fluid Triangulation can be obtained using various software tools For mesh input two files are needed mesh nodes dat and mesh triangles dat The parameters of the mesh are 155 the number of particles on the surface of the immersed object denoted by mesh_nnode and the number of triangular faces in the triangulation denoted by mesh_ntriangle However these parameters are obained automatically from mesh nodes dat and mesh triangles da counting the number of line in respective files The mesh nodes dat thus contains mesh_nnode lines with three real numbers sepa rated by blank space representing three coordinates of the
265. ter setup and would then typically be passed to the function mbtools utils trap_mols mbtools system_generation get_userfixedparts returns the internal list variable userfiredparts which keeps track of all particles that have been fixed in position during the setup This is useful for later releasing particles after warmup routines have been completed mbtools system_generation get_middlebead returns the internal variable middlebead 14 5 2 Available geometries flat fixz bondl arg crystal half pancake shuffle e fizz Fix the vertical positions of all particles The ids of these particles are added to the list of userfiredparts which can later be obtained through a call to mbtools system_generation get_userfixedparts e crystal Sets lipids on a grid instead of randomly 168 half Creates a halfbilayer i e periodic only along one direction Useful to mea sure a line tension pancake Creates a spherical and flat bilayer The diameter of the pancake cannot exceed the box_l shuffle shuffle the topology prior to placing the lipids This is required for a random lipid distribution because otherwise the lipids will be placed on the sphere in the order they appear in the topology Creates a flat bilayer in the XY plane by random placement of lipids sphere c arg initarea arg bondl arg shuffle c 0 0 0 0 0 0 The location of the center of the sphere relative to the center of the box
266. tered through the scripting interface Thus apart from data processing on the fly it can also be used as an efficient correlator for stored data In all cases it produces a matrix of n 2 columns The first two columns are the values of lag times 7 and the number of samples taken for a particular value of 7 The remaining ones are the elements of the n dimensional vector C 7 The uwerr command for computing averages and error estimates of a time series of observables relies on estimates of autocorrelation functions and the respective auto correlation times The correlator provides the same functionality as a by product of computing the correlation function see section 9 2 5 An example of the usage of observables and correlations is provided in the script correlation tcl in the samples directory 9 2 2 Creating a correlation Correlation first has to be defined by saying which observables are to be correlated what should be the correlation operation sampling frequency etc When a correlation is defined its id is returned which is used further to do other operations with the cor relation The correlation can be either updated automatically on the fly without direct user intervention or by an explicit user call for an update Syntax correlation new obs1 id1 obs2 id2 corr_operation operation dt dt tau_max tau_max tau_lin tau_lin compress1 name compress2 name 120 Description Defines a new correlation and returns an integ
267. that the augment argument is optional i e it can be omitted e When an optional argument or a whole command is marked by a superscript label this denotes that the argument can only be used when the corresponding feature see appendix B on page 197 specified in Required features is activated Example 1 constraint wall normal ng ny n dist d type id 2 constraint sphere center Cr Cy Cz radius rad direction direction type id 3 constraint rod center Cs Cy lambda lambda 4 constraint ext_magn_field fr fy fz 2 3 1 ELECTROSTATICS ROTATION DIPOLES 1 Required features CONSTRAINTS 11 2 First steps 2 1 Quick installation If you have installed the requirements see section 1 5 on page 10 in standard locations to compile ESPResSo it is usually enough to execute the following sequence of two steps in the directory where you have unpacked the sources configure make This will compile ESPResSo in a freshly created object path named according to your CPU and operating system As you have not yet specified a configuration a standard version will be built with the most often used features Usually you will want to build another version of ESPResSo with options better suited for your purpose In some cases e g when ESPResSo needs to be compiled for several different platforms or when different versions with different sets of features are required it might be useful to execute the commands not in the source
268. that you read the papers on ELC 6 56 before using it 5 7 7 Dielectric interfaces with the ICC algorithm Syntax iccp3m n_induced_charges convergence convergence_criterion areas areas normals normals sigmas sigmas epsilons epsilons eps_out eps_out relax relazation_parameter max_iterations maz_iterations ext_field ext_field Required features ELECTROSTATICS Description The ICCx algorithm allows to take into account arbitrarily shaped dielectric interfaces This is done by iterating the charge on the particles with the ids 0 to n_induced_particles 1 until the correctly represent the influence of the dielectric discontinuity It relies on a coulomb solver that is already initialized This Coulomb solver can be P3M P3M ELC MMM2D or MMMID As most of the times ICCx will be used with P3M the corre sponding command is called iccp3m Please make sure to read the corresponding articles mainly 7 before using it The particles with ids 0 to n_induced_particles 1 are treated as iterated particles by ICCx The constitute the dielectric interface and should be fixed in space The param eters areas and epsilons are Tcl lists containing one floating point number describing each surface elements area and dielectric constant sigmas allows to take into account a bare charge density thus a surface charge density in absence of any charge induction normals is a Tcl list of Tcl lists with three floating point numbers describing
269. the potentials in the section are therefore singular at zero distance and forces usually become very large for distances below the particle size This is not a prob lem during the simulation as particles will simply avoid overlapping However creating an initial dense random configuration without overlap is often difficult By artificially capping the forces it is possible to simulate a system with overlaps By gradually raising the cap value Fmax possible overlaps become unfavorable and the system equilibrates to a overlap free configuration This command will cap the force to Finar i e for particle distances which would lead to larger forces than Fmax the force remains at Fmax Accordingly the potential is replaced by rFmax Particles placed exactly on top of each other will be subject to a force of magnitude Finax along the first coordinate axis The force capping is switched off by setting Fmax 0 Note that force capping always applies to all Lennard Jones tabulated Morse and Buckingham interactions regardless of the particle types If instead of a force capping value the string individual is given the force capping can be set individually for each interaction The capping radius is in this case not derived from the potential parameters but is given by an additional signal floating point parameter to the interaction 73 6 Setting up the system 6 1 setmd Setting global variables Syntax 1 setmd variable 2
270. the vertices of our lattice which has the spacing a The electric fields E 1 and vector potentials A 1 live on the edges or links and are aligned with them We need also the operator V x It gives the vector B which lives on the faces of the cube or on the plaquettes Fig Figure D 1 Spatial elements of a cell complex In the implementation of the algorithm we assume that particles with masses m and charges q live in the continuum off lattice approach The charges are interpolated on the lattice with grid spacing a using a linear interpolation scheme D 3 Initialization of the algorithm The algorithm as it is implemented only calculates stepwise time updates of the exact field solution Therefore in order to start the simulation for the given random distribution of charges we have to calculate the initial electrostatic field i e the exact solution of the electrostatic problem We find a particular solution of Gauss law as the result of the following recursive procedure see Fig D 2 1 The charge in the plane z Zplane is 1 4 dplane N Y g ri z Zplane D 6 z E 2 N is the number of charges in plane z Zplane Update the z field according to the formula Du 1 dplane E E T epa D 7 205 2 Subtract the charge plane from the each charge on sites of Zplane The charge of the wire Y Ywire Z Zplane S 1 qwire N q ri zi Bn Zplane 9 Yi Ywire D 8 Y i N now mea
271. this is a spanning lipid the first and last elements of this list would typically be head beads e bondtypelist A list of two bondtypeids with the same meaning as explained above for standard lipids Create a lipid which spans across the bilayer 14 5 5 Adding a new molecule type To add a new molecule type you need to define a proceedure which determines how the atoms that make up the molecule should be placed This proc will live directly in the mbtools system_generation namespace Examples can be found in place tcl In order to register your new molecule type to allow placement in any geometry you need to add a call to it in the function mbtools system_generation placemol Make sure that all arguments to your place_mymolecule routine are included in this function call 14 6 Utils Useful utilities routines for various types Includes file management basic geometry and math procedures 173 14 6 1 Setup commands mbtools utils setup_outputdir outputdir paramsfile arg tabdir arg tabnames arg startf arg ntabs arg outputdir Complete path of the directory to be setup At least the parent of the directory must exist e paramfile Name of a file to be copied to the output directory e tabdir Full path name of the directory where forcetables are kept tabnames Complete list of forcetables to be used in the simulation These will be copied to the output directory This routine is designed to setu
272. tic interaction also requires assigning dipole moments to the particles This is done using the part command to set the dipole moment dip e g inter coulomb 1 0 p3m tune accuracy le 4 part O dip 100 part 1 dip 0 O 1 5 8 1 Dipolar P3M Syntax inter magnetic lg p3m Trax mesh cao alpha Required features DIPOLES Description This command activates the P3M method to compute the dipolar interactions between charged particles The different parameters are described in more detail in 10 Teut The real space cutoff as a positive floating point number mesh The number of mesh points as a single positive integer cao The charge assignment order an integer between 0 and 7 alpha The Ewald parameter as a positive floating point number Make sure that you know the relevance of the P3M parameters before using P3M If you are not sure read the following references 13 Note that dipolar P3M does not work with non cubic boxes Tuning dipolar P3M Syntax inter magnetic lg p3m tune tunev2 accuracy accuracy r_cut Tout mesh mesh cao cao alpha a Required features DIPOLES Description Tuning dipolar P3M works exactly as tuning Coulomb P3M Therefore for details on how to tune the algorothm refer to the documentation of Coulomb P3M see section 5 7 1 lor page 62 For the magnetic case the expressions of the error estimate are given in 69 5 8 2 Dipolar Layer Correction DLC Syntax inter magnetic mdlc accur
273. tion that returns a list of floating point numbers of fixed size dimQ can be specified Although its execution might be slow it allows to prototype new observables without a lot of trouble Many existing analysis commands can be made to cooperate with the core analysis that way 9 1 3 Printing an observable Syntax observable id print formatted Description Prints the value of the observable with a given id If the observable refers to the current state of the system its value is updated before printing 9 1 4 Passing an observable to an analysis function Currently the only analysis function which uses the core observables is the correlator section 9 2 9 1 5 Deleting an observable to an analysis function Syntax l observable id delete Description Deletes the observable i e frees the allocated memory and makes the id free for a new observable 9 1 6 Particle specifications You can specify from which particles the observable should be computed in one of the following ways In all cases particle specifications refer to the current state of espresso Any later changes to particles additions deletions changes of types will not be auto matically reflected in the observable 118 e all Requests observable calculation based on all particles in the system e types type_list Restricts observable calculation to a given particle type s The type list is a tcl list of existing particle types e id id list R
274. tive fast linear scaling method for computing induced charges on arbitrary dielectric boundaries J Chem Phys 132 154112 2010 doi 10 1063 1 3376011 S Tyagi A Arnold and C Holm ICMMM2D An accurate method to include planar dielectric interfaces via image charge summation J Chem Phys 127 154723 2007 5 7 3 Sandeep Tyagi Axel Arnold and Christian Holm Electrostatic layer correction with image charges A linear scaling method to treat slab 2d h systems with dielectric interfaces J Chem Phys 129 20 204102 2008 5 7 6 Uli Wolff Monte carlo errors with less errors Comput Phys Commun 156 143 153 2004 222 Index aggregation bending force analysis 97 blockfile Tcl command 127 aggregation blocks bond distances internal first monomer BMHTF interaction 109 bond distances internal first monomer bond lengths 109 center of mass bond lengths chains bond angle interactions end to end distance of a chain bonded interaction type id energies bonded interactions finding holes bonded interactions fsi form factor of a chain box_1 global variable gyration tensor Buckingham interaction hydrodynamic radius of a chain 108 build directory internal distances within a chain 109 local stress tensor 106 cell_grid global variable minimal particle distance cell_size global variable moment of inertia matrix cellsystem Tcl command particle distance center of mass
275. to 0 The total force on a particle can be capped by using the command inter forcecap see section 5 9 5 or on an individual level using the reap variable When Tcap is set and inter forcecap individual has been issued before the maximal force that is generated by this potential is the force at rap By default force capping is off i e the cap radius is set to 0 An optional additional parameter can be used to restrict the interaction from a min imal distance fmin This is an optional parameter set to 0 by default A special case of the Lennard Jones potential is the Weeks Chandler Andersen WCA potential which one obtains by putting the cutoff into the minimum 27 e choosing Teut 260 The WCA potential is purely repulsive and is often used to mimick hard sphere repulsion 5 1 3 Generic Lennard Jones interaction Syntax inter typel type2 lj gen O Teut Cshitt Torr 1 2 bi ba Teaplauto Required features LENNARD_JONES_GENERIC Description This command defines a generalized version of the Lennard Jones interaction see section 5 1 2 between particles of the types typel and type2 The potential is defined by Vat elbiler be 2 Cshift If Tmin Tok lt 1 lt Tout Toft l 0 otherwise 5 2 Note that the prefactor 4 of the standard LJ potential is missing so the normal LJ potential is recovered for b bg 4 ey 12 and ez 6 The total force on a particle can be capped by using the comman
276. ture to the channel given by channelld channelld must be an identifier for an open channel such as the return http www ks uiuc edu Research vmd https github com olenz vtfplugin wiki VTF format 132 value of an invocation of open The output of this command can be used for a stan dalone VSF file or at the beginning of a VTF file that contains a trajectory of a whole simulation Arguments e short verbose Specify whether the output is in a human readable but somewhat longer format verbose or in a more compact form short The default is verbose e radius radii auto Specify the VDW radii of the atoms radii is either auto or a Tcl list describing the radii of the different particle types When the keyword auto is used and a Lennard Jones interaction between two particles of the given type is defined the radius is set to be 4 plus the LJ shift Otherwise the radius 0 5 is substituted The default is auto Example writevsf file radius 0 2 0 1 auto 2 1 0 e typedesc typedesc typedesc is a Tcl list giving additional VTF atom keywords to specify additional VMD characteristics of the atoms of the given type If no description is given for a certain particle type it defaults to name name type type where name is an atom name and type is the type id Example writevsf file typedesc 0 name colloid 1 name pe e ignore_charges this is a temporary workaround for a bug in the VTF reader
277. tzmann fluid while variant 2 turns on the GPU implementation implying that all following LB related commands are executed on the GPU Currently only a subset of the CPU commands are available for the GPU implemen tation For boundary conditions analogous to the CPU implementation the feature LB_BOUNDARIES_GPU has to be activated 12 7 Electrohydrodynamics Syntax setmd mu_E pk uE pE Required features LB LB_ELECTROHYDRODYNAMICS Description If the feature LB_ELECTROHYDRODYNAMICS is activated the non GPU Lattice Boltz mann Code can be used to implicitely model surrounding salt ions in an external electric field by having the charged particles create flow For that to work you need to set the electrophoretic mobility multiplied by the external E field w in all 3 dimensions for your system The three given parameters are float values and should for a meaningful system be less than 1 0 For more information on this method and how it works read the publication 22 153 13 Object in fluid Please cite 11 BIBTEX key cimrak in file doc ug citations bib if you use the object in fluid implementation described below Simulations using ESPResSo work mostly with objects molecules atoms polymers colloids crystals that are physicaly composed of points linked together with bonds These objects are like skeletons without inner or outer volume The idea is to use ESPResSo for objects that do have i
278. uces to the textbook result min 1 exp 1 T gt 2 1 T Ua Up kg 7 2 However T can also be chosen to be any other parameter for example the Bjerrum length i e the the strength of the electrostatic interaction In this case 6 T is a constant but the energy Uc T of a configuration C depends on T and one needs the full expression 7 1 ESPResSo always uses this expression In practice one does not swap configurations but temperatures simply because ex changing temperatures requires much less communication than exchanging the properties of all particles Th ESPResSo implementation of parallel tempering repeatedly propagates all config urations C and tries to swap neighboring temperatures After the first propagation the routine attempts to swap temperatures T and To T3 and T4 and so on After the second propagation swaps are attempted between temperatures To and 73 T4 and T5 and so on For the propagation parallel tempering relies on a user routine typically one will simply propagate the configuration by a few 100 MD time steps Details on usage and an example The parallel tempering code has to be loaded explicitely by source scripts parallel_ tempering tcl from the Espresso directory To make use of the parallel tempering tool one needs to implement three methods the propagation the energy calculation and an initialization routine for a configuration A typical initialization routine will look
279. uch of the user time while trying to understand what the code is supposed to do Hand in hand with the extensibility and readability of the code comes the flexibility of the whole program On the one hand it is provided by the generalized functionality of its parts avoiding highly specialized functions An example can be the implementation of the Generic Lennard Jones potential described in section where the user can change all available parameters Where possible default values are avoided providing the user with the possibility of choice ESPResSo cannot be aware whether your particles are representing atoms or billiard balls so it cannot check if the chosen parameters make sense and it is the user s responsibility to make sure they do On the other hand flexibility of ESPResSo stems from the employment of Tcl at the steering level Apart from the ability to modify the simulation and system parameters at runtime many simple tasks which are not computationally critical can be implemented at this level without even touching the C kernel For example simple problem specific analysis routines can be implemented in this way and made interact with the simulation core Another example of the program s flexibility is the possibility to integrate system setup simulation and analysis in one single control script ESPResSo provides commands to create particles and set up interactions between them Capping of forces helps prevent system blow up whe
280. ume View change_volume Inew x y z xyz velocities Umax start pid count N Ses ef SRL S e s fp 2441313 invalidate_system parallel_tempering main rounds N swap swap perform perform init init values T connect master port port load jnode resrate Nreset info info parallel_tempering set_shareddata data metadynamics metadynamics metadynamics metadynamics metadynamics metadynamics metadynamics metadynamics metadynamics Required features METADYNAMICS analyze mindist type_list_a type_list_b analyze distto pid analyze distto Y 2 analyze nbhood pid r_catch analyze nbhood x y z reatch analyze distribution part_type_list_a part_type_list_b rmin rmaz rbins log_flag int_flag analyze radial_density_map zbins ybins xrange yrange axisofrotation centerofrotation beadtypelist thetabins analyze modes2d analyze get_lipid_orients set off set distance pid pidg dmin dmax bheight bwidth Joouna dbins set relative_z pid pidg min Zmax Oheight width bound bins print_stat current_coord print_stat coord_values print_stat profile print_stat force load_stat profile_list force_list ele 2 St sf sf y OSS S analyze lipid_orient_order analyze bilayer_set analyze bilayer_density_profile analyze cell_ accuracy numberofinteractions analyze get_folded_positions molecule shift zx y z analyze Vkappa reset read set Vi Vig avk analyze
281. ument tmar defines the maximum value of t for which G r t is calculated If it is omitted or set to zero maximum possible value is used If the particles perform a random walk i e a normal diffusion process G r t r is a Gaussian distribution for all times Deviations of this behavior hint on another diffusion process or on the fact that your system has not reached the diffusive regime In this case it is also very questionable to calculate a diffusion constant from the mean square displacement via the Stokes Einstein relation Output format The output corresponds to the blockfile format see section 10 2 on page 127 msd msd 0 msd 1 vanhove G 0 0 G 1 0 G 0 1 G 1 1 The G r t are normalized such that the integral over space always yields 1 8 1 14 Center of mass Syntax analyze centermass part ype Description Returns the center of mass of particles of the given type 8 1 15 Moment of inertia matrix Syntax 1 analyze momentofinertiamatrix typeid 2 analyze find_principal_axis typeid 101 Description Variant 1 returns the moment of inertia matrix for particles of given type typeid The output is a list of all the elements of the 3x3 matrix Variant 2 returns the eigenvalues and eigenvectors of the matrix 8 1 16 Gyration tensor Syntax analyze gyration_tensor typeid Description Analyze the gyration tensor of particles of a given type typed or of all part
282. ut Toff inter typel type lj cos2 O Tog w Required features LJCOS LJCOS2 inter typel type smooth step o n kg 02 Teut Required features SMOOTH_STEP inter type type2 bmhtf nacl A B C Do Te Required features BMHTF_NACL inter type type morse Q Tmin Teut Required features MORSE inter typel type buckingham A B C D Teut Taiscont shift Required features BUCKINGHAM inter typel type soft sphere a Nn Teut Toffset Required features SOFT_SPHERE inter typel type hat Fmax Te Required features HAT inter typel type2 hertzian o Required features HERTZIAN inter typel type gaussian O Teut Required features GAUSSIAN inter typel type lj angle O Teut bla blo b2a b2 reap 20 Oz k Required features LJ_ANGLE inter typel type gay berne o Up Teutor k1 k2 y y Required features ROTATION GAY_BERNE inter bondid fene K Armax ro inter bondid harmonic K R reut inter bondid subt_lj reserved R inter bondid rigid_bond constrained_bond_distance positional tolerance velocity_tolerance inter bondid tabulated bond filename inter bondid tabulated angle filename inter bondid tabulated dihedral filename inter bondid virtual_bond inter bondid stretching_force L p ks inter bondid bending_force 0 kp inter bondid area_force_local Sipo kal inter bondid area_force_global 9 kag inter bondid volume_force V ky E si sf l l Si Si Sf l l SISPSPey Es w alakale 1 TE
283. values of p x If p x contains values larger than 1 default value of max the maximum or any number larger than that has to be given maz This routine basically takes the function p x and places it into a rectangular area 0 1 0 max Then it uses to random numbers to specify a point in this area and checks wether it resides in the area under p x Attention Since this is written in tcl it is probably not the fastest way to do this vec_random len returns a random vector of length len uniform distribution on a sphere This is done by chosing 3 uniformly distributed random numbers 1 1 If the length of the resulting vector is lt 1 0 the vector is taken and normalized to the desired length otherwise the procedure is repeated until succes On aver age the procedure needs 5 739 random numbers per vector This is probably not the most efficient way but it works Ask your favorit mathematician for a proof phivec_random lt v gt lt phi gt len return a random vector at angle phi with v and length len 141 e PARTICLE OPERATIONS Operations involving particle positions The parameters pi can either denote the particle identity then the particle position is extracted with the The part command command or the particle position directly When the optional box parameter for minimum image conventions is omited the functions use the the setmd box_1 command bond_vec lt p1 gt lt p2 gt Calculate bond vector pointi
284. warmup may not be possible due to physical restrictions ESPResSo cannot detect these mistakes and it is your responsibility to find simulation procedure suitable to your specific problem 19 3 Getting compiling and running ESPResSo This chapter will describe how to get compile and run the ESPResSo software ESPResSo releases are available as source code packages from the ESPResSo home pagd This is where new users should get the code The code within release packages is tested and known to run on a number of platforms Alternatively people that want to use the newest features of ESPResSo or that want to start contributing to the software can instead obtain the current development code via the version control system software git from ESPResSo s project page at the Savannah GNU server This code might be not as well tested and documented as the release code it is recommended to use this code only if you have already gained some experience in using ESPResSo Unlike most other software no binary distributions of ESPResSo are available and the software is usually not installed globally for all users Instead users of ESPResSo should compile the software themselves The reason for this is that it is possible to activate and deactivate various features before compiling the code Some of these features are not compatible with each other and some of the features have a profound impact on the performance of the code Therefore it is not poss
285. when using the development sources doxygen Creates the Doxygen code documentation in the doc doxygen subdirectory tutorials Creates the ESPResSo tutorials in the doc tutorials subdirectory doc Creates all documentation in the doc subdirectory only when using the develop ment sources A number of options are available when calling make The most interesting option is probably j num_jobs which can be used for parallel compilation on computers that have more than one CPU or core num_jobs specifies the maximal number of jobs that will be run Setting num_jobs to the number of available processors speeds up the compilation process significantly 3 3 Running ESPResSo When ESPResSo is found in your path it can be run via Espresso tcl_script Largs 23 When ESPResSo is called without any arguments it is started in the interactive mode where new commands can be entered on the command line When the name of a tcl_ script is given the script is executed Any further arguments are passed to the script If you want to run ESPResSo in parallel using MPI the actual invocation depends on your MPI implementation In many cases e g OpenMPI the command will be mpiexec n n_nodes Espresso tcl_script args where n_nodes denotes the number of MPI nodes to be used However note that de pending on your MPI installation MPI jobs can only be run in a queueing system so that ESPResSo will not run from the command line Also old
286. writes the interactions in a standardized format to channel 10 2 4 Writing the random number generator states Syntax 1 blockfile channel write random 2 blockfile channel write bit_random 3 blockfile channel write seed 4 blockfile channel write bitseed Description Variants 1 and 2 write the full information on the current states of the respecitive random number generators see sections 11 2 1 on page 144 and 11 2 2 on page 145 on any node to channel Using this information it is possible to recover the exact states of the generators Variants 3 and 4 write only the seed s which were used to initialize the random number generators Note that this information is not sufficient to restore the full state of a random number generator because the internal state might contain more information 10 2 5 Writing all stored configurations Syntax blockfile channel write configs Description This command writes all configurations currently stored for off line analysis see sec tion 8 3 on page 112 to channel 129 10 2 6 Writing arbitrary blocks Syntax 1 blockfile channel write start tag 2 blockfile channel write end 3 blockfile channel write tag arg Description channel has to be a Tcl channel Variant 1 starts a block and gives it the title tag variant 2 ends the block Between two calls to the command arbitrary data can be written to the channel When variant 3 is used the function
287. y nz radius rad length length type id constraint rod center Cr Cy lambda lambda 7 1 7 8 constraint plate height h sigma sigma 9 constraint ext_magn_field fr fy fz 2 3 10 constraint plane cell z y z type id 11 constraint mindist_position x y z 1 Required features CONSTRAINTS ELECTROSTATICS ROTATION DIPOLES Description The constraint command offers a variety of surfaces that can be defined to interact with desired particles Variants 1 to 6 create interactions via a non bonded interaction potential where the distance between the two particles is replaced by the distance of the center of the particle to the surface The constraints are identified like a particle via its type for the non bonded interaction After a type is defined for each constraint one has to define the interaction of all different particle types with the constraint using the inter command In variants 1 to 5 constraints are able to be penetrated if flag is set to 1 Otherwise when the penetrable option is ignored or flag is set to 0 the 36 constraint cannot be violated i e no particle can go through the constraint surface In variants 1 to 4 it is also possible to specify a flag indicating if the constraints should be reflecting The flags can equal 1 or 2 The flag 1 corresponds to a reflection process where the normal component of the velocity is reflected and the tangential component remains unchanged If the flag is 2
288. y the partial slip boundary condition 6B Onvlrp Ulea where vj denotes the tangential component of the velocity and Op its spatial derivative normal to the surface both evaluated at the position rg of the so called hydrodynamic boundary This boundary condition is characterized by two effective parameters namely i the slip length dg and ii the hydrodynamic boundary rg Within the approach of the tunable slip boundary interactions it is possible to tune the slip length systematically from full slip to no slip A coordinate dependent Langevin equation describes a viscous layer in the vicinity of the channel walls which exerts an additional friction on the fluid particles T is the temperature yz the friction coefficient and Teut is the cut off radius of this layer t is the timestep of the integration scheme With vr vy and v it is possible to give the layer a reference velocity to create a Plane Couette Flow Make sure that the cutoff radius rey is larger than the cutoff radius of the constraint Lennard Jones interactions Otherwise there is no possibility that the particles feel the viscous layer This method was tested for Dissipative Particle Dynamics but it is intended for meso scopic simulation methods in general Note that to use tunable slip boundary interac tions you have to apply two plane cell constraints with Lennard Jones in addition to the tunable slip interaction Make sure that the cutoff radius reut is larger than the
289. y the particle triples p1 pa and p3 and p2 p3 and pa Together with appropriate Lennard Jones interactions this potential can mimic a large number of atomic torsion potentials If you enable the feature OLD_DIHEDRAL then the old less general form of the potential is used V K 1 p cos n 5 32 where p is rather a phase factor and can only take values p 1 60 5 7 Coulomb interaction Syntax 1 inter coulomb 0 0 2 inter coulomb 3 inter coulomb parameters Description These commands allow to set up the calculation of the Coulomb interaction The Coulomb or electrostatic interaction is defined as follows For a pair of particles at distance r with charges q and q2 the interaction is given by UC r Ink ar EE 5 33 where lg e2 4rekgT denotes the Bjerrum length which measures the strength of the electrostatic interaction As a special case when the internal variable temperature is set to zero the value of bjerrum length you enter is treated as lgkgT rather than lg This occurs when the thermostat is switched off and ESPResSo performs an NVE integration see also Section 6 2 Computing electrostatic interactions is computationally very expensive ESPResSo features some state of the art algorithms to deal with these interactions as efficiently as possible but almost all of them require some knowledge to use them properly Unedu cated use can result in completely unphysical simulations V
290. ze and the centre of mass position is calculated for each block separately com_velocity particle_specifications blocked size Velocity of the centre of mass If blocked size is specified the particles are subdi vided into blocks of size size and the centre of mass velocity is calculated for each block separately 116 com_force particle_specifications blocked size Total force on the specified particles If blocked size is specified the particles are subdivided into blocks of size size and the total force is calculated for each block separately stress_tensor The stress tensor It only works with all particles It is returned as a 9 dimensional array Ora Oxy Oxz yx Oyy Oyz Ozx Ozy Ozz stress_tensor_acf_obs The observable for computation of the Stress tensor autocorrelation function Sim ilarly to the stress tensor it only works with all particles It is returned as a 6 dimensional array Ory Oyz Ozz Ore Oyy Oxa Ozz Oyy Ozz where oj are the components of the stress tensor particle_currents particle_specifications Electric currents due to individual particles For a particle i j qu At where At is the simulation time step Required feature ELECTROSTATICS currents particle_specifications Electric currents averaged over all particles j gt qjv7 At where At is the simulation time step Required feature ELECTROSTATICS dipole_moment particle_specifications The dipole m

ESPResSo User's Guide

Contents

Download Pdf Manuals

Related Search

Related Contents