Downloads - CompChemMPI
Contents
1. 2802 picturename strrep titlename picturename strrep picturename title picturename FontSize 18 xe lab kl tb p Fonts ize 218 ylabel rho kg m 3 FontSize 18 set gca FontSize 18 Lexend a b c Location Best NA 29 Code 40 analysingenergy m 2 Matlab savename strrep fnames 3 K 1 J name 2 xvg png savename2 strcat path Figures density savename saveas f savename2 png fclose a fprintf 30s 3s 10 2f 3s 10 2f 3s 10 2f 3s n titlename amp mean dens11 2000 10000 mean dens12 2000 10000 mean dens 3 2000 10000 XV tmp mean dens 1 2000 10000 mean dens 2 2000 10000 mean dens13 2000 10000 1 meanmean mean tmp errormean std tmp A prontf 42 90s423s210 2149s 10 1 f 23s4n totlenameys E y meanmean errormean AN end 30 20 21 22 23 24 25 26 27 28 2D maps e g orientation analysis Andrey Frolovs scripts for plotting the 2D maps are given here The script PlotOrientationIL m calls the function CF_orient_2D stored in CF_orient_2D m Code 41 CF_orient_2D m Matlab function CF orient 2D x phi cdatati ZCREATEFIGURE CDATA1 2 1 image cdata fs 30 figurel figure XVisual 70x27 Tru2. 23 24 25 26 27 28 Code 13 topol top ASCII file define _FF_OPLS define _FF_OPLSAA defaults nbfunc comb rule gen pairs fudgeLJ fudgeQQ 1 3 yes 0 5 0 5 333 LOAD ATOM TYPES include path top EMI_AtTy_lopes itp include path top BMI_AtTy_lopes itp include path top OMI_AtTy_lopes itp include path top Cl_AtTy itp include path top BF4_AtTy itp include path top TFSI_AtTy_lopes itp 333 LOAD OPLS FF include localgromacspath share gromacs top oplsaa ff ffnonbonded itp include localgromacspath share gromacs top oplsaa ff ffbonded itp 333 LOAD MOLECULES itp include path top EMI_lopes itp include path top BMI_lopes itp include path top OMI_lopes itp include path top Cl itp include path top BF4 itp include path top TFSI lopes itp system Name Neat BMI BF4 molecules BMI 200 BF4 200 1 2 3 Charges In the case of organic molecules the OPLS AA forcefield has proved to give reasonable results in many cases But it might be a good idea to play with the charges With Gaussian03 3 charges can be calculated based on quantum mechanics An input file com for calculating charges of a molecule would look like Code 14 Head of Gaussian input com ASCII file ANprocshared 4 Mem 1GB Chk scna_HF_6 31Gd chk tp hf 6 31g d nosymm geom connectivity pop chelpg with the addition of particle positions and bond information 11 20 21 2
3. Code 30 g rdf non interactive Bash script bin bash for k in 0 1 2 do g rdf f traj k xtc s NPT k tpr n rdf index ndx o rdf k Cation Anion xvg rdf mol com b 2000 lt lt EOF Cation Anion EOF g rdf f traj k xtc s NPT k tpr n rdf index ndx o rdf k _Cation_SOL xvg rdf mol com b 2000 lt lt EOF Head SOL EOF done 22 3 2 3 Distribution function in cylindrical geometry In the case of cylindrical geometry the distribution function around a CNT can be calculated using the flagg xy if the nanotube is oriented in z direction Code 31 rdp sh Command line g rdf bin bin b begin e end f run xtc s run tpr n index ndx o fname cn cn_ fname rdf mol com xy The resulting rdfs will not go to one in infinite space but rather to a constant value different from one This is due to the solvent excluded volume of the CNT and therefore the rdf has to be renormalized 3 2 4 Distribution function in slab geometry Take a look on the Gromacs manual for Interface related items e g g_order g_density g potential g traj In addition to those tools it is possible to use the Gromacs tool g_rdf with the flagg surf taken the carbon atoms of the graphene wall as reference The tool will caclulate the distance between any atom center of mass with respect to the atoms of the surface This is not completly correct as the length of the hypotenuse is not identical with th
4. PBS_O_WORKDIR Shift to the directory we submitted the job from cd PBS_O_WORKDIR Load the CASTEP module module add ze gromacs run msd sh with run_msd sh Code 24 run msd sh Bash script bin bash for i in ls d EMI do cd i echo Cation g_msd f traj xtc s NVT tpr n rdf_index ndx o msd_Cation xvg b 4000 echo Anion g_msd f traj xtc s NVT tpr n rdf index ndx o msd Anion xvg b 4000 cd done 18 2 4 3 German supercomputer JUROPA For running simulations on the German supercomputer JUROPA the following start script might be useful submitted using msub runscript_JUROPA Code 25 runscript_J UROPA pbs Bash script bin bash MSUB l nodes i ppn 8 MSUB l walltime 24 00 00 MSUB v tpt 1 module load gromacs 4 5 5 module load mk1 10 2 5 035 Prepare the Gromacs run grompp f NVT mdp c packmol gro p topol top n indez ndz o NVT Run the MD on 8 cores mpiexec np 8 exports GMXLIB GROMACS_ROOT bin mdrun_d deffnm NVT maxh 24 Further details can be obtained from the JUROPA webpage with quick introductions http www fz juelich de ias jsc EN Expertise Supercomputers JUROPA UserInfo QuickIntroduction html 2 5 Comments e When running simulations with periodic molecules it is very important to include in the mdp file the line periodic molecules yes otherwise e g graphene sheets will form a ball e For slab geo
5. R Cheeseman Montgomery Jr J A T Vreven K N Kudin J C Burant J M Millam S S Iyengar J Tomasi V Barone B Mennucci M Cossi G Scalmani N Rega G A Petersson H Nakatsuji M Hada M Ehara K Toyota R Fukuda J Hasegawa M Ishida T Naka jima Y Honda O Kitao H Nakai M Klene X Li J E Knox H P Hratchian J B Cross V Bakken C Adamo J Jaramillo R Gomperts R E Stratmann O Yazyev A J Austin R Cammi C Pomelli J W Ochterski P Y Ayala K Morokuma G A Voth P Salvador J J Dannenberg V G Zakrzewski S Dapprich A D Daniels M C Strain O Farkas D K Malick A D Rabuck K Raghavachari J B Foresman J V Or tiz Q Cui A G Baboul S Clifford J Cioslowski B B Stefanov G Liu A Liashenko P Piskorz I Komaromi R L Martin D J Fox T Keith M A Al Laham C Y Peng A Nanayakkara M Challacombe P M W Gill B Johnson W Chen M W Wong C Gonzalez and J A Pople Gaussian 03 Gaussian Inc Wallingford CT 2004 W Humphrey A Dalke and K Schulten VMD Visual molecular dynamics Journal of Molecular Graphics 14 1 33 38 February 1996 35
6. o 150 H 0 108 150 0 0 140 o and should be extended by the last two lines of the following piece Code 10 atomname2type n2t new ASCII file C opls_145 0 12 12 011 3 C 0 150 C 0 150 H 0 108 C opls_145 0 12 12 011 3 C 0 133 C 0 150 0 0 140 C opls_145 0 0 12 011 3 C 0 140 C 0 140 C 0 140 C opls_145 0 0 12 011 2 C 0 140 C 0 140 The positions of CNT atoms are restrained to the initial values by harmonic potential with 1000 kJ mol nm force constant in each direction For restraining the positions of the carbon atoms the Gromacs tool genrestr can be used It generates an itp file that has to be included in the general topology file of the nanocarbon molecule 1 2 2 Organic molecules In the upper section it was mentioned that the inital configuration of organic molecules were taken from databases or self prepared with the help of SCHR DINGER Maestro software 9 0 211 1 Again this program is used to generate apropriate OPLS AA parameter See Code 3 on page 7 how to start Maestro After building the molecule alternatively reading in a coordinate file downloaded from a database the next step is to create a cms file with system coordinates and force field applied This is done by using the menu Applications Desmond System Builder In TA nice tutorial as given here http chembytes wikidot com gromacs wiki 20 21 22 23 24 the window that pops up the solvent model shou
7. Gromacs tools One reason for the wide acceptance and usage of Gromacs is the enormous number of analysis tools that are provided by the developer and which can be found at http www gromacs org Documentation Gromacs_Utilities They are all preinstalled and can be started just by using the command line as it is possible for grompp or mdrun 3 2 1 Density and density profile The first look that should be taken after a simulation finished sucessfully and also if not is how temperature energy volume etc evolved over time This can be done straightfoward using the tool g_energy Code 27 g_energy Bash script bin bash for Cation in EMI BMI OMI do for Anion in Cl BF4 TFSI do cd Bulk_neat Cation _ Anion for k in 0 1 2 do echo 7 8 9 10 11 16 17 g_energy f ener k edr s NPT k tpr o Cation _ Anion _ k xvg done cd done done The resulting xvg files can be read in by numerous programs e g Gnuplot Matlab Grace and analyzed further e g plotted statistical analysis of mean std etc See section 4 on page 27 20 Density profiles can be obtained by using the Gromacs tool g_density which also allows to calculate charge densities but unfortunately not number densities of a center of mass Therefore it might be more convenient to use g_rdf but with using different flaggs See the subsection below for the usage 3 2 2 Radial distribution function and coordination number Again straigh
8. called anneal Then it greps coordinate files out of these trr files using the Gromacs tool trjconv The resulting frames are stored with the given name name _rep gro plus arunning number e g name rep0 gro name _rep1 gro name rep2 gro etc Code 19 grepreplicas sh Bash script bin bash for itrr in ls anneal x xtc do name echo itrr sed s xtc g echo 0 n q triconv s name tpr f name xtc o name _rep gro sep b 120 dt 50 done 2 3 3 System heating After the preparation of one coordinate file it is possible to heat the system using the annealing algorithm Code 20 annealing ASCII file annealing single annealing_npoints 2 annealing_time 0 1000 annealing_temp 1500 350 15 In this example two annealing points are set the system starts at time 0 ps with a temperature of 1500 K and within the next 1000 ps the system is smoothly cooled down to the simulation temperature of 350 K In a system with fast diffusing particles this proceadure is very useful but the highter the viscosity of a system the less happens during the heating and therefore resulting systems cannot be taken as independent 2 4 Computing facilities 2 4 1 Local systems For starting Gromacs on local systems with more than one core and if gromacs is installed as a parallel version one should always use the flagg mdrun nt 1 for using only one thread or how many threa
9. ei br a OS s ssy B y 6 Starting Maestro Command line ee ee 7 genbox Command line oy A is A 3 7 packmolimpurity inp ASCII file cir dia ra 8 sed Command ais Es AE AD rr e ee a 8 Starting Packmol Command line 4 0 a pra 2 8 hops Command line ea e a ee a ea ae ea 9 atomname2type n2t ASCII file kree oc ake oe w n Oe a vin ee 9 ahommame2ivpen2t new ASCII file e dy A eee b Po 9 ASCE fle iTe aga A Oak tw R N 10 BE ip Ad A AS 10 topoltopr ASGIT ley ls Moi o Rat b de TAN 11 Head of Gaussian input com ASCII file 11 INE Imdp i ASCH file e aca ae bd eo a 12 NPT mdpA2 ASC file x oo oy de ee Ag 13 grompp Command line te gat ye ok hk ak pe A Tk Don 13 mart OMIA iba to Gh si a disp e e b o ee XKS 13 grepreplicas sh Bash script 2444 222 44 book a bee ba be dae 203 15 annealing O i n kot s Sent od beer m es son e RO l 15 mdrun nt 1 Command line cs S en ee See ew YA ee x 16 runseript parallel pbs Bash script wa 6a ka eh oes hee ee 17 runscript serial pbs Bash seript n hades ds ete eee Shr ae ee sak 18 run msd sh Bash seript ih ios ey an eh atte ds Gave rien Gulla haven e eka Bohs Deeg 18 runscript JUROPA pbs Bash script tka sil won m CO a 19 Start VMD Command bo ep abla As die Be ale Peed akan ls 20 psenerov dha dng be A anan BE aoe OS Dee
10. further analysis a Gat ee Sea eee Running replicas ac ot asas a G s b ni z 2 3 1 2 3 2 2 3 3 Independent initial configurations by different random seeds Grabbing frames from a trajectory System MOAN fon te ec Seen g doh wen at dot eG DG Aye Pa Compimnetaci tes z ya ae ale ce ee din he o Geers Bayete 2 4 1 2 4 2 2 4 3 Local systems wes te Mwe al nerds AN aa Se Dees British supercomputer HECTOR tea tb German supercomputer JUROPA 2 0 AC OMMIMMEMGS or is oe e ea GA ee Mee ae AA ety nc hee Oy deen Sp Ge ub System analysis 3 1 Vis alanalysis on a e ake pata a Se ee A a aaa de poke oe ap ee ey sw 9 2 GrOMACSItOOIS ma e aes wR oe Bo eos ees WAS eae ee 3 3 3 2 1 3 2 2 3 2 3 3 2 4 3 2 5 3 2 6 3 2 7 3 2 8 Density and density profile y tr Ye ae be GA ave de ANA e Am Radial distribution function and coordination number Distribution function in cylindrical geometry Distribution function in slab geometry Order parameter for alkyl chains br A a Order parameter for alkyl chains as a function of box length Head stacking How to analyse the orientation of aromatic rings 2D number density map xx x d SA Self written O ya qa aea sal ef s de eee 3 3 1 3 3 2 3 3 3 3 3 4 Potential of mean force and free energy Orientation anal
11. semiisotropic tau p 1 0 1 0 compressibility 4 5e 5 0 0 1le 5 ref_p 1 0 1 0 bar GENERATE VELOCITIES FOR STARTUP RUN gen vel yes gen temp 298 1 gen seed 473529 OPTIONS FOR BONDS constraints hbonds constraint_algorithm lincs unconstrained_start no shake_tol 0 00001 lincs_order 4 lincs_warnangle 30 morse no lincs_iter 2 2 2 Running Gromacs In general it is necessary to produce a run input file for Gromacs first This is done by the tool grompp which also helps detecting numerous input errors Code 17 grompp Command line grompp f 1 NPT mdp c steep gro p topol top o NPT Then the molecular dynamics run can be started with mdrun Code 18 mdrun Command line mdrun deffnm NPT maxh mdrun deffnm NPT maxh mdrun deffnm NPT maxh mdrun deffnm NPT maxh path2 path3 path4 cpi NPT cpt append multi 4 1 1 1 1 cpi NPT cpt append multidir pathi 13 The flagg deffnm saves all files under the given name but with proper extension NPT log NPT trr NPT xtc and so on This is especially useful when starting in the same folder several simulations like steep NPT and NVT The flagg maxh t stops Gromacs after 0 99 t hours while writing a checkpointfile for the last step The second line shows how to continue the sim ulation using the checkpoint file with cpi NPT cpt and appending the new output to the old files using appe
12. the simulations are liable to crash These limits are due to network comunication 14 2 2 3 Data storage for further analysis In case of water simulations a rule of thumb tells to sample coordinates each 0 3ps In case of room temperature ionic liquids 1 ps seems already sufficient and might be due to the slow dynamics increased even more if disc space is an issue For analysis of velocity autocorrelation function and ionic conductivity the velocities of the system should be sampled each 40 fs 0 04 ps according to Maginn for 4ns be cautious this fast results in several GB of disk space needed per simulation 2 3 Running replicas There are three possibile ways to get a set of independet system configurations replicas for improving analysis quality and covering the whole thermodynamic ensemble 2 3 1 Independent initial configurations by different random seeds Preparation of independent initial configurations by using different seeds with Packmol Seed numbers that have been tested so far are 191917 171719 151517 191317 111719 Any large primary number should do the job as well 2 3 2 Grabbing frames from a trajectory For simulating replicas it might be useful to take the first sucessful production run and only grep a few configurations out of it While assigning random velocities and or additional heating in many cases the configurations can be taken as independant The following script takes the xtc files from a directory
13. 1 1 2 Combining them all to full systems under study If the coordination files xyz pdb or gro of all single compounds of the system under study have been created the preparation of the complete simulation box is straightforward It can be done by either using the gromacs tools like genbox or by using the free software PACKMOL 2 Using genbox would look like Code 4 genbox Command line genbox cp graphene sheets gro cs C6mimBF4 gro maxsol 1240 with having a predefined configuration of two graphene sheets graphene_sheets gro and a file with a preequilibrated simulationbox of neat solvent C6mimBF4 gro Due to the packing mechanism placing the whole bulk in the empty space and removing all overlapping molecules this method demands time computational effort and has no safety that it really works It is NOT recommended The most convinient way to produce a working simulation box is to use PACKMOL 2 which uses a mathematic geometric filling algorithm rather than an overlapping test of hard spheres Packmol creates an initial point for molecular dynamics simulations by packing molecules in defined regions of space The packing guarantees that short range repulsive interactions do not disrupt the simulations The great variety of types of spatial constraints that can be attributed to the molecules or atoms within the molecules makes it easy to create ordered systems such as lamellar spherical or tubular lipid layer
14. 2 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 2 Running simulations with Gromacs 2 1 Molecular dynamics parameter file mdp After creating preparing the coordinate and topology files the only missing files are the mdp files which define the simulation parameter like integrator annealing proceadure tem perature etc A sample file is given here Code 15 NPT mdp 1 ASCII file RUN CONTROL PARAMETERS integrator md tinit 000 dt 0 002 nsteps 500000 comm mode Linear nstcomm 1 energy_grps EMI TFSI OUTPUT CONTROL OPTIONS nstxout 1000 nstvout 1000 nstfout 0 nstlog 1000 nstenergy 50 nstxtcout 500 xtc_precision 1000 xtc_grps NEIGHBORSEARCHING PARAMETERS nstlist 10 ns_type grid pbc xyz periodic_molecules yes rlist 1 3 OPTIONS FOR ELECTROSTATICS AND VDW coulombtype PME rcoulomb 143 rcoulomb switch 1 0 vdw_type Shift rvdw 14 0 fourierspacing 0 12 pme_order 4 ewald_rtol 1 05 ewald_geometry 3dc optimize_fft lt o n 12 20 21 22 23 24 25 26 27 28 29 Code 16 NPT mdp 2 ASCII file OPTIONS FOR WEAK COUPLING ALGORITHMS tcoupl v rescale tc grps System tau_t 1 0 ref_t 298 1 Pcoupl Berendsen Pcoupltype isotropic tau_p 0 compressibility 4 5e 5 1le 5 ref_p 1 0 bar Pcoupltype
15. 5 code 15 is shifted to the end of file This decreases readability a lot Codes that are naturally longer than one page have to be splitted by hand into two different listing environments minipage requires a complete new environment definition with the possible loss of caption information and unexpected behaviour newpage Introduce white space and needs to be rechecked after every change of the document Codes that are naturally longer than one page have to be splitted by hand into two different listing environments e Changed basicstyle footnotesize to basicstyle ttfamily in the header of the ETEX document option of lstset to avoid overlapping of symbols in code examples e g and M in the Gaussian input header 18 April 2012 e In subsubsection German supercomputer JUROPA JUROPA runscript added e In section Running Gromacs flaggs added multi multidir e In section Plotting fitting and statistical analysis subsection General comments on figures added e In subsection Plotting data subsubsection Grace Xmgrace added e In subsection Gromacs tools subsubsection Order parameter for alkyl chains added In subsection Gromacs tools subsubsection Order parameter for alkyl chains as a function of box length added e In subsection Gromacs tools subsubsection 2D number density map added e In subsection Gromacs tools subsubsection Head stacki
16. Manual Computational Physical Chemistry SCRIPTS and TOOLS for and INFORMATION on Running and Analysing Molecular Dynamics Simulations Kathleen Kirchner Physics and Life Sciences Nanoscience Division Department of Physics Strathclyde University G4 ONG Glasgow U K April 14 2012 Short summary This document should provide the reader with detailed information on starting and analysing molecular dynamics simulations with GROMACS It is a collection of scripts and experiences that have been made by Andrey Frolov and Kathleen Kirchner under the supervision and with the help of Prof Dr Maxim Fedorov starting from the middle of 2008 begin of Andrey Frolov s PhD till now Contents 1 System preparation Tl WOM A AA TON a Y nc im Boe ee we os Be a a ee ae ee ee 4 1 2 1 1 1 1 1 2 Generation of first coordinates ciones a tet Combining them all to full systems under study Force fields 2400 wey min kin aaa 77 1 2 1 1 2 2 1 2 3 Carbon nanoparticles s i Aa b dei si b ea ay e Re eS aai Organic molecules 2 2 5a Dt an Dt eo R RK oe eae A Charges oe Bx orn os O ek A oe ec 2 Running simulations with Gromacs 2 1 Molecular dynamics parameter file mdp 2 2 Te AE AN Ss IN pee bad so da a A 2 3 2 4 2 2 1 22 2 2 2 3 Energy minimization and equilibration Production runs ooa a s s e ee De Data storage for
17. be used to estimate the free energy of the process of dividing one water molecule from a full solvated ion by calculating the depth of the first minimum of the ion water PMF 3 3 2 Orientation analysis The Gromacs tool g_sorient analyzes solvent orientation around solutes It calculates two angles between the vector from one or more reference positions to the first atom of each solvent molecule the angle between a vector spanned by 3 defined atoms A1 A2 A3 and the vector spanned by another atom B1 plus the first atom of the predefined atoms A1 Therefore the atoms of the solvent to use have to be specified using an index file On the preparation of index files see Code 29 on page 22 Modified version of g sorient for cylindrical geometry Orientation distributions are calculated with the g_sorient which has to be modified to be able to handle cylindrical symmetry of CNTs and to be able to resolve the orientation probability density as a function of distance The usage is the following Code 37 g sorient Command line g_sorient f runi xtc s runi tpr n index ndx o sori_OMI_COM xvg ro sord MI C M xvg xy com b 2000 e 32000 rmax 3 0 nat2 37 ati2i 13 ati22 9 ati2cen 1 rbin 0 01 cbin 0 1 The distributions of cos 04 for rmin lt r lt rmaz are calculated with respect to the center of mass of the solvent molecules com and along the z axis xy The following list describes all additional options implem
18. ds one wishes to use Otherwise the program will occupy everything Code 21 mdrun nt 1 Command line grompp f O_STEEP mdp c packmol k gro p topol top o steep k maxwarn 1 mdrun nt 1 deffnm steep k grompp f 1 NPT mdp c steep k gro p topol top o NPT k mdrun nt 1 deffnm NPT k amp 16 20 21 22 23 24 25 2 4 2 British supercomputer HECTOR For running simulations on the British supercomputer HECToR the following start script might be useful submitted using qsub runscript parallel pbs Code 22 runscript_parallel pbs Bash script bin bash login PBS N NAME PBS q parallel PBS l mppwidth 24 PBS l mppnppn 24 PBS A z01 fedo PBS V PBS l walltime 01 00 00 l cput 00 05 00 echo PBS_O_WORKDIR Shift to the directory we submitted the job from cd PBS_O_WORKDIR module add xe gromacs Get the number of MPI tasks and tasks per node export NPROC qstat f PBS_JOBID grep mppwidth awk print 3 export NTASK qstat f PBS JOBID grep mppnppn awk print 31 tpr NPT MAXH 1 aprun n NPROC N NTASK mdrun_mpi maxh MAXH deffnm tpr dlb auto cpi state cpt append 17 Analysis e g g rdf can be done in serial using the following submission script Code 23 runseript serial pbs Bash script bin bash login PBS N NAME PBS q serial PBS A z01 fedo PBS V PBS l cput 01 00 00 echo
19. e The following script reads in all energy xvg files obtained through echo 7 8 9 10 11 16 17 g energy f ener edr s NPT tpr o energy xvg in the current directory and plots the column eight vals 8 density versus column one valsf1 time This is done for three replicas on one goal which have seperate file names _a xvg _b xvg and c xvg In addition the mean value of a certain part hopefully a plateau is calculated and printed in a format directly usable in 28 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 Code 39 analysingenergy m 1 Matlab clear all close all fnames dir xvg numfids length fnames fprintf 30s 10s 3s n RTIL amp Impurity amp Nrho kg m 3 fprintf 30s 3s 10s 3s 10s 3s 10s 3s n JO DOT KOU AN for K 1 numfids 3 for J 1 3 count 3 K 1 J a fopen fnames count name rt vals textscan a tit eye he Le LE NE NE s Headerlines 8 CommentStyle Q time J vals 1 dens J vals 8 end f figure visible otf plot time 1 dens 1 time 2 dens 2 time 3 dens 3 LineWidth 3 tmp strrep fnames 3x K 1 41 name _c xvg tmp strrep imp _ tmp strrep tmp BF4 BF 4 tmp strrep tmp 1M no neat titlename strrep tmp SOL
20. e exact distance between a flat wall and the atom but in most cases the induced error will be less than 3 assuming a distance between carbon atoms of the wall of 1 4 Aand a distance between atom and wall of at least 5 A Code 32 dpz sh Command line bin bash g rdf bin bin b begin e end f run xtc s run tpr n index ndx o fn cn cn_ fn rdf rdftype surf mol nopbc EOF 1 p2 EOF 23 3 2 5 Order parameter for alkyl chains g_order allows the calculation of the order parameter for alkyl chains angle between z axis and vector spanned by 3 distinct carbon atoms Preparation of a proper index file is required the resulting index file should contain only entries of the carbon atoms that belong to the alkyl chain In the case of imidazolium cations also the first Nitrogen is added to the list Code 33 index_order ndx ASCII file No 2 21 40 59 78 97 116 135 154 173 192 211 230 249 268 1 64 1 13 32 51 70 89 108 127 146 165 184 203 222 241 260 279 C2 9 28 47 66 85 104 123 142 161 180 199 218 237 256 275 3 2 6 Order parameter for alkyl chains as a function of box length Use g_order s1 to calculate the order parameter as a function of box length The general proceadure is explained in the section above 3 2 7 Head stacking How to analyse the orientation of aromatic rings g sgangle z calculates the angle between z axis and vector spanned by 3 defined atoms Unf
21. eColor depth 24 RGB mask Oxff0000 Oxff00 0x001f De ZZfigurel1l fiqure XVisual Oxz27 TrueColor depth 24 RGB mask Ozff0000 Ozff00 Ox00ff set figurel Position 1300 300 800 6001 1 axes Parent figure1 LineMidth 2 Layer top FontSize fs ONT Teles 0 5 0 Ob A so XTick 0 5 1 0 1 5 1 7 2 0 2 5 gt XMinorTick on gt YMinorTick on TicklL ngih 0 02 006 si CLim 0 1 4 Uncomment the following line to preserve the X limits of the ares box axesi on hold axes1 7all 4 Create image image x phi cdatal Parent axes1 CDataMapping scaled xlim axes1 10 6 1 71 2 Uncomment the following line to preserve the Y limits of the azes ylim axes1 1 1134 caxis 10 1 51 xlabel r nm FontSize fs ylabel cos phi FontSize fs colorbar peer axes1 FontSize fs 31 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 Al 42 43 Code 42 PlotOrientationIL m Matlab path path pud cur dir pwa x 0 00 0 01 2 99 phi 0 95 0 1 0 95 list dir 4 2M ACN EMI TFSI 2 ACN TFSI 2M TFSI TFSI F3 list_sub_dir 1 q 057 PQQ y 2g05 r mols_list GAT ea TESI rs for d 1i len
22. eadure with the resulting configuration several times Andrey Frolov figured out while working on energy minimization that it is not sufficient to run Gromacs in single precision to obtain reasonable results Convergence of the energy is in most cases only reached when using Gromacs in double precision this is done by recompiling the source code For system equilibration NPT simulations should be performed In cases of low viscous flu ids it takes a few 10 ps until the density reaches a constant value In case of room temperature ionic liquids it might be necessary to equilibrate at least 1 ns 2 2 2 Production runs A production run means a molecular dynamics simulation run that is resulting in enough data to sample the ensemble space correctly and that allows statistical reliable analysis of data The length of a production run depends on the system size the more molecules are included in the simulation box the highter is the probability so sample all configurations within a given time frame The length also depends on the viscosity of a system or how fast particles are moving equal to how fast they forget about theirs past The systems under study in our group demanded for simulation lenghts between 30 ns and 100 ns One option to ensure sufficient data is the usage of replica runs as explained in Section 2 3 on page 15 SGromacs shows reasonable parallelization if not less than 1000 atoms per core are used Below 500 atoms per core
23. ented by Andrey Frolov 25 g sorient c flaggs C nat1 FALSE etINT amp nat1 Number of atoms in the reference molecules nat2 FALSE etINT amp nat2 Number of atoms in the molecules to calculate orientation atil1 FALSE etINT amp ati11 Vector origin Atom index on reference molecule ati12 FALSE etINT amp ati12 Vector end Atom index on reference molecule ati21 FALSE etINT amp ati21 Vector origin Atom index on molecule to calculate orientation ati22 FALSE etINT amp ati22 Vector end Atom index on molecule to calculate orientation ati2cen FALSE etINT amp ati2cen Vector end Atom index on molecule to calculate orientation The following description of the analysis routine was taken from Andrey Frolov s PhD thesis We calculated the average number of particles at a certain distance r around CNT in the following way NAr 0 01 nm r po 2nr Ar 0 where po is the bulk density of the corresponding particles me is the RDP of the corresponding particles NAr 0 011m 7 is the average number of particles at a certain distance r in the cylindrical volume segment with the difference between radii of smaller and larger cylinders of Ar To use the modified Gromacs script it is necessary to recompile the complete source code That means the follwing steps have to be done e Compile the Gro
24. gth list dir for sdel length list sub dir for n i length mols list if strcmp list dir d 2M ACN BMI TFSI mols_list 1i BMI elseif strcmp list dirfdi 2M A N UMI TFS I mols_list 1i OMI else mols_list 1 EMI end fn sori_ mols listin _COM dat dfn cur_dir list_dirfd list_sub_dir sd J arr load dfn c reshape arr 3 20 3001 CF_orient_2D x phi c xlim 0 6 1 7 caxis 10 1 51 ofn cur_dir list_dirfd list_sub_dirf sd mols _listin p g 1 ofn ofn syek ks eval print dpng r200 ofn 1 close all end end end exit a Ja fn 32 4 3 Fitting data 4 3 1 Matlab Fitting in Matlab can be done using splines An amazing code can be downloaded from Matlabs file exchange http www mathworks com matlabcentral fileexchange 13812 It is called splinefit and smoothes even the noisiest data set 4 4 Statistical analysis ToDo 4 4 1 Matlab 4 42 R 5 How does your system look like when As force fields are only model approaches of real systems they may not cover the correct physical behaviour in all imaginable circumstances The formation of crystals is one problem in simulating diluted systems 5 1 it crystalizes The simulation normally do not crash but nevertheless do not show a behaviour that has been observed in experiments To assure the formation of crystals there are two straightfor
25. he common part of the simulation tools developed at CAMd ASE provides Python modules for manip ulating atoms analyzing simulations visualization etc Code 1 graphenesheet py Python from ase import from ase structure import graphene_nanoribbon gnri graphene nanoribbon 3 4 type armchair sheet True cell gnri get_cell print cell posx 11 0 0 posy cell 1 1 posz cell 2 2 io urite graphene nanoribbon 3 4 armchair pdb gnr1 The resulting pdb file can be transformed rotated etc using editconf Code 2 editconf Command line editconf s graphene nanoribbon 3 4 armchair pdb translate 0 rotate 0 0 0 o confout trans rot gro N methyl 2 pyrrolidone NMP acetonitrile AN as vvell as tetraethylammonium TEA tetrabutylammonium TBA and imidazole based ionic liquids like 1 ethyl butyl octyl 3 methylimidazolium EMIm BMIm OMIm with the anions Cl tetrafluoroborate BF4 and bis trifluoromethylsulfonyl imide TFSI have been modelled using the OPLS AA force field http turin nss udel edu research tubegenonline html 2http www ccp14 ac uk ccp web mirrors jcrystal products wincnt 3https wiki fysik dtu dk ase epydoc ase structure module html The inital configuration were taken from databases or self prepared with the help of SCHRO DINGER Maestro software 1 Code 3 Starting Maestro Command line SCHRODINGER maestro SGL
26. is 20 gedi Command line cscs cored oats oa eee Bek ee ee he 21 make ndx non interactive Bash script wade OE Ee EG 22 g rdf non interactive Bash seript wo eo oa eed ket eh ad ees 22 TAPE Commend line sss o a e e a 23 d zsl Command line ct aia di kk S soul 23 index_order ndx ASCII file ya A AI A AA 24 g sgangle z workflow A ork x Scan bok Sok dwet bed eo Bebe s eS 24 PSOA OLE WOR EHOW cas ie tn bata ty BR ape Da an he de iat be te gee the RAP 24 g densmap Command line he ae a RA 25 g sorient Command line Gan ee de Ee a e ee 25 Compiling Gromacs Command line 26 analysingenergy m 1 Matlab o A eS eg 29 analysingenergy m 2 Matlab o o eee 30 Matlab Ste oa ood ke as x yx ee ea 31 PlotOrientationIL m Malabia der ee S wee a ee 32 Recent changes 14 April 2012 e Several newline added to prevent codes from being splitted over two pages NOTE This should only be a temporary solution Possibilities to prevent latex from splitting the code on two pages Istlisting environment options float ht The main problem is that the code is embedded in a float environment but with a fixed starting point This fixed starting point can be deleted by the option begin lstlisting float ht This will hold the code on one page but also shift the code somewhere e g code 6 appears before code
27. ld be none all other entries can be set as default No ions should be added Now press start The resulting desmond_setup out cms file can be then transfered to Gromacs files by using the script Maestro2gmx py To do Usage of the script Experience showed that it is very convenient to store all topologies itp files in an own directory and separate atomtypes with mass charge and non bonded interaction parameter Lennard Jones form Code 11 BF4_AtTy itp ASCII file atomtypes type mass charge ptype sigma epsilon B 10 811 0 8276 A 3 5814e 01 3 9748e 01 F 18 998 0 4569 A 3 1181e 01 2 5104e 01 and bonded interation parmeter in case of inconsistency the upper charges in _AtTy itp are used Code 12 BF4 itp ASCII file moleculetype name nrexcl BF4 3 atoms nr type resnr residu atom cgnr charge mass 1 B 1 BF4 B 1 0 8276 10 811 2 F 1 BF4 F 1 0 4569 18 998 3 F 1 BF4 F 1 0 4569 18 998 4 F 1 BF4 F 1 0 4569 18 998 5 F 1 BF4 F 1 0 4569 18 998 bonds ai aj funct co ci 1 2 1 0 137 284512 000 1 3 1 0 137 284512 000 1 4 1 0 137 284512 000 1 5 1 0 137 284512 000 angles ai aj ak funct 0 1 3 1 2 1 110 611 502 080 4 1 2 1 110 611 502 080 4 1 3 1 110 611 502 080 5 1 2 1 110 611 502 080 5 1 3 1 110 611 502 080 5 1 4 1 110 611 502 080 The simulation directory should contain then a top file where all necessary or even more itp files are included 10 20 21 22
28. macs code on a local folder Code 38 Compiling Gromacs Command line 1 configure prefix your gmx folder 2 make e Copy the gmx sorient c to the your gmx folder src tools folder and remove the old files Removing old files Command line irm g sorient g_sorient o e Recompile Gromacs Recompiling Gromacs Command line 1 cd your gmx folder 2 make 26 Now the executable your gmx folder src tools g_sorient is new Orientation probability density in 2D maps The modified version of g_sorient produces a solvent orientation map sord_ xvg with the name provided by the flagg g_sorient ro It can be read in and visualized using Matlab and the image object 3 3 3 Volumes of solute cavities ToDo The volumes of cavities of molecules can be calculated with the help of Gaussian03 software 3 In this example the geometries of the species are optimized at the B3LYP 6 31g d p level of theory In the quantum mechanics calculations the Self Consistent Isodensity Polarizable Continuum Model SCI PCM is used to model the acetonitrile solvent To do add script 3 3 4 Residence time ToDo Old script written in FORTRAN by Andrey Frolov Will take some time to get the key points 4 Plotting fitting and statistical analysis 4 1 General comments on figures Always check title axis labels and units before submitting any figure to anyone In general most plotting tools allow st
29. metries levald geometry 3dc should be added but this line has to be avoided when doing bulk simulations The simulations may either crash or show weired physical behaviour such as formation on vacuum bubbles or continous increasing and decreasing of volume breathing e On some computers simulations are crashing due to some error using PME dynamic load balancing Therefore by using the option mdrun dlb no the dynamic load balancing can be switched off and systems may run more stable 19 3 System analysis 3 1 Visual analysis The program of choice for looking at trajectories is Visual Molecular Dynamics VMD 4 It is useful to read in a trajectory track special particles atoms molecules residues etc and prepare nice pictures of the molecular systems They provide also analysis tools for example of radial distribution functions or mean square deviation from an input structure but if possible it is preferable to use GROMACs tools or self written scripts A Gromacs trajectory can be read in via the terminal by using the following command this is much more preferable than reading in a trajectory using the gui for the trajectory is not displayed and therefore reading in is much faster Code 26 Start VMD Command line vmd gro test gro xtc test xtc The box boundaries can be shown using Extensions Tk Console and typing in pbc box To do Add some useful lines for particle selection 3 2
30. nd Two usefull flaggs multi nand multidir pathi path2 path3 path4 with the same goal are introduced in line 3 and 4 of the above starting lines for mdrun Gromacs allows to summarize several independend simulations into one job allowing even for small systems the usage of super computers 5 In case of multi nn tpr files have to be stored in one folder being numbered as follows NAMEO tpr NAME1 tpr NAME2 tpr NAME3 tpr The flagg multidir 1 path2 path3 path4 was introduced in Gromacs version 4 5 4 unfortuantely an entry in the general help mdrun h is still missing With multidir tpr files can be stored in differend folders path but should all have the same name NAME tpr It depends on the personell preferences which of these two options to use 2 2 1 Energy minimization and equilibration The first step after system preparation should be a basic energy minimization to remove high forces that would lead to a system explosion In Gromacs this is done by changing the integrator to steep in the mdp file using integrator steep The number of steps can range between 1000 and several 10000 If the initial configuration was created with Packmol and the initial density was reasonable only a few minimization steps are sufficient to relax the system and start molecular dynamics simulations If a bigger number of minimization steps is necessary it also might be better to do energy minimizations of 1000 steps and redo the proc
31. ng How to analyse the ori entation of aromatic rings added http www mathworks com matlabcentral fileexchange 13812and http web cecs pdx edu gerry MATLAB plotting loadingPlotData htmladded as useful Matlab scripts 1 System preparation For all substances under study initial coordinates xyz and gromacs topology files itp and AtTy itp are stored in the subfolder ForceFields In addition a Collection of GRO MACS topologies for small organic molecules developed and maintained by David van der Spoel Sweden and Carl Caleman Germany can be found at http virtualchemistry org BENCH and GROMACS user contribution for topologies at http www gromacs org Downloads User_contributions Molecule_topologies 1 1 Configurations 1 1 1 Generation of first coordinates Initial configurations of single wall carbon nanotubes CNT can be generated using the online tool TubeGen It demands the input of chirality and number of replications and generates for example a pdb file with the resulting structure Nanocarbon onions can be modelled as a collection of carbon fullerenes of different size e g C720 C320 C60 The coordinates of these substructures can be obtained through databases like the Fullerene Library by M Yoshida or special tools like the Nanotube Modeler Graphene layers can be again prepared in several ways One version is the use of ase structure a tool of the Atomic Simulation Environment ASE that is t
32. orage of the figures in many formats most handy is eps as it can be read in by To convert figures to a specific format there are in general three options 1 Use an image editor like Gimp Very slow in case of many figures 2 A simple convert figure jpg figure eps in the shell should do the trick on unix systems 3 More advanced is the usage of Inkscape inkscape file figure jpg export eps figure e Numerous options like flipping are explained here http tavmjong free fr INKSCAPE MANUAL htm1 CommandLine html 4 2 Plotting data 4 2 1 Grace Xmgrace The software Xmgrace is distributed with the Gromacs software or can be downloaded using the synaptic package manager keyword Grace The simple command Starting Xmgrace Command line Xmgrace NAME xvg 27 allows to visualize any xvg output file The possibilities to change the representation are limited therefore plotting with Matlab or Gnuplot might be more handy for high quality images The advantage of Xmgrace lies in its fast usage and the possibility to make labels and units given in the header of xvg human readable To make an eps file use File Print setup and choose as device EPS This only sets up the printing do get the eps file printed use File Print 4 2 2 Matlab A handy tool for plotting data with Matlab can be downloaded at http web cecs pdx edu gerry MATLAB plotting loadingPlotData html xvg e g density vs tim
33. ortunately only one molecule per calculation is allowed This will introduce the following procedure for an analysis of all cations Code 34 g sgangle z workflow While Browse through all cations do Make index file for cation do Run g_sgangle z do Sum up histogram of angle distribution Print normalized histogram g sgangle calculates the angle between 2 vector spanned by 6 defined atoms The proce dure gets somewhat more complicated for the analysis of all cations Code 35 g sgangle workflow While Browse through all cations do Make index file for cation While Browse through all other cations do Make index file for other cation do Run g_sgangle do Sum up histogram of angle distribution Print normalized histogram 24 3 2 8 2D number density map An interesting tool for visualizing number density distribution is g densmap Check carefully the range that is taken into account when analysing interfaces Code 36 g densmap Command line g densmap f NVT xtc s NVT tpr aver z xmin 1 65 xmax 2 15 b 4000 e 10000 convert densmap xpm densmap eps 3 3 Self written scripts 3 3 1 Potential of mean force and free energy The calculation of the potential of mean force PMF can easily be done by using any math ematical programming language e g Matlab FORTRAN taking the initial obtained radial distribution functions and calculating PM F r kg T In g r The PMF can
34. s The user must provide only the coordinates of one molecule of each type the number of molecules of each type and the spatial constraints that each type of molecule must satisfy http www gromacs org Documentation How tos Mixed_Solvents Shttp www ime unicamp br martinez packmol 7 2 A working packmol input would look like Code 5 packmol impurity inp ASCII file tolerance 2 0 filetype pdb output packmol pdb seed seednum add_box_sides structure top All_itp Cation pdb number 200 inside box 0 0 0 50 0 50 0 50 0 end structure structure top All_itp Anion pdb number 200 inside box 0 0 0 50 0 50 0 50 0 end structure structure top All_itp Impurity pdb number numberofmol inside box 0 0 0 50 0 50 0 50 0 end structure A second script was prepared for changing the name of Cation Anion and Impurity as well as the number of impurity molecules the box size and the seed number for generating independent replicas of the system using the stream editor sed Seed numbers that have been tested so far are 191917 171719 151517 191317 111719 Any large primary number should do the job as well Code 6 sed Command line sed s 50 0 Box g s Impurity Imp s numberofmol numImp s Cation Cation s Anion Anion s seednun seed packmol impurity inp gt packmol tmp inp The final steps are starting PACKMOL and after a successful r
35. tforward is the calculation of the radial distribution function g r with g_rdf The general usage is Code 28 g rdf Command line g rdf f traj xtc s NPT tpr en rdf_index ndx o rdf Cation Anion xvg rdf mol com b 2000 with an interactive input of the groups that should be used By default these groups are the whole system and one group for every molecule type In this case it is phsical reasonable to use the flagg rdf mol_com to calculate the rdfs between the center of masses and not some cummulative value between all atoms which is done by VMD s tool so be careful 21 If the task is to analyse rdfs between certain atoms it is necessary to use new groups by creating a customized index file which is done by the Gromacs tool make_ndx Code 29 make ndx non interactive Bash script bin bash make ndx f confoutO gro o rdf_index ndx lt lt EOF keep 0 ri 1 200 name 1 Cation ri 201 400 name 2 Anion a OW a HWi a HW2 name 3 SOL a Nia1 a Nibi a C4al a C4b1 a C4c1 a Nia3 a Nib3 a C4a3 a C4b3 a C4c3 a Nia4 a Nib4 a C4a4 a C4b4 a C4c4 name 4 Head a Cibi a Cid3 a Cid4 name 5 Tail q EOF With the given index file g rdf can be used again this time without interactive usage but rather within a shell script Notice the usage of group names instead of numbers which makes the usage of g_rdf less prone to error and reusable if the system structure changes
36. un the transformation of packmols output pdb file to the gromacs input file gro with the help of editconff Code 7 Starting Packmol Command line HOME Programs packmol packmol lt packmol_tmp inp editconf f packmol pdb o packmol gro Shttp manual gromacs org current online editconf html 1 2 Force fields 1 2 1 Carbon nanoparticles The generation of the force field for carbon nanoparticles is described in detail in Andrey Frol ovs tutorial on simulating carbon nanotubes CNTs Tutorials Tutorial simulate CNT Non bonded interaction parameters for the nanotube nanoonion graphene carbon atoms cor respond to the benzene OPLS AA all atom optimized molecular potential for liquid simu lation carbon opls 145 in Gromacs notation This was done by using the Gromacs tool x2top Code 8 x2top Command line x2top f CNT gro o CNT top name C60 kb 392459 2 kt 527 184 pbc The flagg name renames the molecule default is ICE In the case of periodic carbon nanotubes and graphene sheets it is very important to add the flagg pbc to assure that bonds between all particles are recognized Be aware the file ffoplsaa n2t Gromacs version 3 x or atomname2type n2t Gromacs version 4 x needs an aditional line for recognizing all carbon bond The following two lines are already included Code 9 atomname2type n2t ASCII file C opls_145 0 12 12 011 3 C 0 150 C C opls_145 0 12 12 011 3 C 0 133 C
37. ward ways that can be performed in an early stage of work and not that expensive in terms of computational time First visualising the trajectory should be done If only the ions are shown one can possibly observe the formation of pairs chains and later agglomerates Second the piecewise calculation of the radial distribution function between ions can be done Andrey Frolov performed the simulation of 0 15 M Nal in NMP for 100ns He cutted the trajectory into pieces of 10ns and caclulated for every piece the radial distribution functions between Iodide atoms The resulting figure showed a continious increase of the peak size indicating the stepwise formation of an agglomerate of particles Be aware that the effect gets more and more pronounced when simulation time increases So simulations should run at least for 60 ns 6 Useful script lines ToDo 6 1 Shell 6 1 1 Improvements for the bashrc alias 18 15 color auto alias l ls color auto export PATH scripts nanotube PATH 33 6 1 2 Bash 6 1 3 sed cat awk 6 2 Matlab 6 3 Python 6 4 Fortran 34 References 1 2 Maestro 9 1 User Manual May 2010 L Mart nez R Andrade E G Birgin and J M Mart nez Packmol A package for build ing initial configurations for molecular dynamics simulations Journal of Computational Chemistry 30 13 2157 2164 October 2009 M J Frisch G W Trucks H B Schlegel G E Scuseria M A Robb J
38. ysis urea s ki ae las A fe Ot Volumes of solute cavities ToDo Residence time ToDo 258000 les oe ee ler k S 12 12 13 14 14 15 15 15 15 15 16 16 17 19 19 4 Plotting fitting and statistical analysis 4 1 General comments on figures oa bce Pa SSO LAS CLAY HESS 4 2 Plotting data seii ie oa o Poa e A S KE 50 5 lit 4 2 1 Grace Xmgrace Fok p s A dam k ye b oe a ee yon 25007 Matlab s A A 0 A MATTEO data e chy eG oda AR eb bsy a b s Se wi ee be 200 0 ab Eb ez b ozun o o ds a A A S 4 4 Statistical analysis ToDo ao ea Ee OR daa R dan at ie 4 LA Matlaba sey a Me y a x RE A a Mo eki YR riy Be IN dena BPE eq S Oe Dan dae af v 6 Useful script lines ToDo Oy Hell eismas p Bake Se Pate wore M h m d y eet OE 6 1 1 Improvements for the bashre 20 paternidad ra Bl Bash c 3 s n 0 bb b 109 GS educa o o ES y E yoo A A A s ees ee BM 0597 PUNT De agate As Ete Mee ba bb b da See wiles b ip Ora se IA S de SA 27 27 27 27 28 33 33 33 33 33 Codes CONDO WN FER H Co Co Co Co Co Co Qo CO A Qo No BHR m RFR HE O Ot He Co OO 1 Ob Ot HS Co CD O OO Ob Ot He Co graphenesheet py Python wi a wa wate an dik v OT deya WOR we 6 editconf Command nel
Download Pdf Manuals
Related Search