Tutorial-1.md

@@ -233,8 +233,58 @@ execute cpptraj with
 
 The resulting file 02_Heat_re.mdcrd can now again be viewed with vmd (see figure)
 
+## STEP3: Production run 
 
+At this stage, we will let the system equilibrate at 300 K for 600 ps to get a realistic starting configuration for the QM/MM calculations. This simulation uses constant pressure, thus the volume of the octahedral box will adapt. The input is again given below:
 
+![step3__pressure](uploads/fc30233bda00c8aaaa3073ecaf761c62/step3__pressure.png)
+
+Start the computation with:
+
+`$AMBERHOME/bin/sander -O -i 03_Prod.in -o 03_Prod.out -p BPH_solv.top -c 02_Heat.rst -r 03_Prod.rst -x 03_Prod.mdcrd -inf 03_Prod.mdinfo`
+
+and analyse the output of the dynamics using the amber utility:
+
+`mkdir 03_Prod; cd 03_Prod; process_mdout.perl ../03_Prod.out`
+
+going back one directory and plotting both runs into one plot with xmgrace
+
+`cd ..
+xmgrace 02_Heat/summary.TEMP 02_Heat/summary.EKTOT 02_Heat/summary.EPTOT \ 02_Heat/summary.ETOT 02_Heat/summary.PRES 03_Prod/summary.TEMP \ 03_Prod/summary.EKTOT 03_Prod/summary.EPTOT 03_Prod/summary.ETOT \ 03_Prod/summary.PRES`
+
+
+yields two areas (zoomed in Figure 5) with the values of interest. In the lower section (-32000 to -20000) we find potential and total energies, the upper area (-1000 to 800) contains pressure and temperature info. Below 5000 is the kinetic energy (not plotted here). The heating process can be followed by the raise in all values between 0 and 20 ps. After ca. 50 ps simulation time, the system oscillates around the mean values of energy, temperature and pressure values with no significant drop or increase, so that we can safely assume that it has equilibrated by the end of the simulation. The summary.VOLUME values further indicate that the initial volume has dropped and also oscillates around the mean value. You can also watch the results of the simulation in VMD by modifying the “re-image.ptraj” script and running cpptraj, following the same steps as described in STEP 2. For the QM/MM simulation, we will take a point close to the end of the simulation which is near to the average value of the potential energy.1 For my simulation, step number 2976 (with EPTOT=-25157.7397) is in this range (search for a corresponding point in your calculation). For further processing,  we will extract a spherical droplet with waters in a distance of 20 Å from the BPH center with cpptraj. 
+
+1 Normally you would rather sample over several snapshots and average the resulting values, but this is a functionality tutorial. More details on sampling is given in the main manual.
+
+![figura5](uploads/f5941e57fb6ec137bf3e4e42ba7e6a18/figura5.png)
+
+Figure 5: Results of Amber MD calculations plotted with xmgrace.
+
+The input file for taking the snapshot with cpptraj (snap.trajin) is documented below. We will directly take the mdcrd file so we have to re-image first.
+
+`trajin 03_Prod.mdcrd 2976 2976 1`  # read from mdcrd frames 2976 to 2976 (1 frame)
+`center :1` 					# put BPH in the center
+`image familiar` 				# re-image 
+`trajout snapshot.rst rest` 			# write the coordinates of this frame
+`go`					  # proceed
+
+To delete all water molecules with a distance > 20 A from the center (we will use C7 of BPH here) we need another script called “strip.in”
+
+
+`trajin snapshot.rst rest`			# read the snapshot
+`reference snapshot.rst rest` 	# use it as reference (necessary for strip command)
+`strip @7>:20.0`				# strip all waters further than 20A around atom C7 
+`trajout strip.pdb pdb`			# write pdb output
+`go`						# proceed
+
+execute with
+
+`cpptraj -p BPH_solv.top -i snap.trajin`
+`cpptraj -p BPH_solv.top -i strip.trajin`
+
+
+The resulting PDB file strip.pdb (see Fig. 6) is the basis for our QM/MM setup