3.1 Large scale simulation using the parallel I/O
Files for this tutorial (tutorial-3.1.tar.gz)
Contents
- 3.1.1 Making restart for minimization
- 3.1.2 Minimization with restraints on the protein
- 3.1.3 Making restart for heat-up
- 3.1.4 Heat-up with restraints on the protein
- 3.1.5 Making restart for equilibration 1
- 3.1.6 Equilibration with restraints on the macromolecules
- 3.1.7 Making restart for equilibration 2
- 3.1.8 Equilibration without restraints
- 3.1.9 Making restart for production
- 3.1.10 Production
- 3.1.11 Convert
Let us see the contents in the downloaded tutorial files by using the following UNIX commands. The input files to be used in this tutorial are already prepared in the tutorial-3.1/toppar directory. In this directory, we have villin_crowding.pdb
(PDB file), villin_crowding.psf
(protein structure file), par_all36_prot.prm
(CHARMM parameter file), top_all36_prot.rtf
(CHARMM topology file), and toppar_water_ions.str
(STR file).
# make a new directory to do this tutorial $ mkdir tutorial-3.1 $ cd tutorial-3.1 # copy and decompress the downloaded file $ cp ~/Download/tutorial-3.1.tar.gz ./ $ tar -xzf tutorial-3.1.tar.gz $ cd tutorial-3.1 $ ls 1_prst_min 2_min 3_prst_heat 4_heat 5_prst_eq1 6_eq1 7_prst_eq2 8_eq2 9_prst_prod 10_prod 11_convert toppar
Note: The parallel I/O scheme are available in SPDYN only.
In this tutorial, we explain an example of the parallel I/O MD simulation for macromolecular crowding system, which contains multiple proteins (villin headpiece).
The procedures of the simulation are identical with conventional ones, except that we are dealing
with the parallel I/O using the same number of input files and the number of MPI processes.
In all cases, we use 8 MPI processes.
3.1.1 Making restart for minimization
To make use of the parallel I/O, the first step to do before minimization is making the parallel restart files, each of which will be assigned to each MPI process for minimization run. prst_setup included in genesis-install-directory/bin is used to make parallel restart files.
The control file for prst_setup is very similar to that for minimization except for a few changes. First, in [OUTPUT] section, instead of output restart file, parallel restart file for minimization should be specified (the brackets corresponds the indexes of the MPI ranks) . Second, you do not need [MINIMIZE] section (you can just leave the same information written in the control file of the minimization). Third, you should specify the number of domains in each direction or the total number of domains (which is identical to the total number of MPI processes). Below, Input (prst_min.inp) to generate multiple restart files for minimization is shown. It is difficult to distinguish the necessary keywords from unnecessary ones. So, we wrote all the keywords identically to control file for minimization. It is necessary to change in [OUTPUT] and [BOUNDARY] sections only. In [BOUNDARY] section, you need to set the number of sub-domains in each direction using domain_x, domain_y, and domain_z, or total number of sub-domains using domain_xyz. In this example, we set the number of sub-domains in each direction:
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str # CHARMM stream file psffile = ../toppar/villin_crowding.psf # CHARMM protein structure file pdbfile = ../toppar/villin_crowding.pdb # PDB file reffile = ../toppar/villin_crowding.pdb # PSF file [OUTPUT] rstfile = ./setup().rst # GENESIS restart file cachepath = ./cache keep_cache = NO [ENERGY] pairlistdist = 13.5 # pairlist distance [SELECTION] group1 = an:N|an:CA|an:C|an:O # selection group 1 [RESTRAINTS] nfunctions = 1 # number of restraint function 1 function1 = POSI # [DIST,POSI,RMSD] [MASS] kind of function 1 constant1 = 10.0 # force constant(x,y,z) for function 1 select_index1 = 1 # restraint groups in function 1 [BOUNDARY] box_size_x = 150 # box size (x) in [PBC] box_size_y = 150 # box size (y) in [PBC] box_size_z = 150 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
Another noticeable change is cachepath in [OUTPUT] section. It is a path to a directory where intermediate files are written before generating final restart files. If you do not specify cachepath, all the intermediate files will be saved in memory only, and it takes more time. After executing these commands, you generate 8 restart files named setup0.rst ~ setup7.rst. These files will be used instead of PDB and PSF files for minimization.
3.1.2 Minimization with restraints on the protein
After making multiple restart files for minimization, you can minimize with 8 MPI processes, as the number of restart files is 8. Each MPI process read one corresponding restart file (i.e., MPI rank 0 reads setup0.rst). The control input file of minimization does not need PSF and PDB files.
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str rstfile = ../1_prst_min/setup().rst [OUTPUT] rstfile = run().rst # restart file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [MINIMIZE] method = SD nsteps = 1000 rstout_period = 1000 eneout_period = 10 nbupdate_period = 10 force_scale_init = 0.00005 force_scale_max = 0.0001 [SELECTION] group1 = an:N|an:CA|an:C|an:O [RESTRAINTS] nfunctions = 1 # number of functions function1 = POSI # restraint function type constant1 = 10.0 # force constant select_index1 = 1 # restraint group [BOUNDARY] box_size_x = 150 # box size (x) in [PBC] box_size_y = 150 # box size (y) in [PBC] box_size_z = 150 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.3 Making restart for heat-up
In this process, [CONSTRAINTS] section, which does not exist in the previous control input file is newly added. Unlike the restart files without the parallel I/O, the potential function information is included in the restart files (that’s why PSF file is not necessary for the parallel I/O). So any change of potential functions requires a consistent restart files to be regenerated. The procedure of regenerating restart files is very similar to the case of minimization. The control input file to generate restart files is as follows:
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str psffile = ../toppar/villin_crowding.psf # CHARMM protein structure file pdbfile = ../toppar/villin_crowding.pdb # PDB file reffile = ../toppar/villin_crowding.pdb rstfile = ../2_min/run().rst [OUTPUT] rstfile = setup().rst # DCD trajectory file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NVT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat [SELECTION] group1 = an:N|an:CA|an:C|an:O [RESTRAINTS] nfunctions = 1 # number of functions function1 = POSI # restraint function type constant1 = 10.0 # force constant select_index1 = 1 # restraint group [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for x-axis domain_z = 2 # # of domain for x-axis
3.1.4 Heat-up with restraints on the protein
Except for the [INPUT] section, the control input for heat-up is the same as those without parallel I/O.
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str rstfile = ../3_prst_heat/setup().rst [OUTPUT] dcdfile = run().dcd # DCD trajectory file rstfile = run().rst # restart file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] integrator = LEAP # [LEAP,VVER] nsteps = 5000 # number of MD steps timestep = 0.002 # time step (ps) eneout_period = 1 # energy output period crdout_period = 500 # coordinates output period rstout_period = 5000 # restart output period annealing = YES # simulated annealing anneal_period = 50 # annealing period dtemperature = 3.4 # temperature change at annealing (K) [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NVT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat temperature = 0.1 # initial temperature (K) [SELECTION] group1 = an:N|an:CA|an:C|an:O [RESTRAINTS] nfunctions = 1 # number of functions function1 = POSI # restraint function type constant1 = 10.0 # force constant select_index1 = 1 # restraint group [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.5 Making restart for equilibration 1
For equilibration, we take two-step process in this tutorial. First, we equilibrate only the solvent with positional restraints on the protein backbones. Second, the restraints will be removed, and protein structures are relaxed in the solvent. Because [ENSEMBLE] and [RESTRAINTS] sections are changed, we should regenerate the restart files for each process using prst_setup. The control input to generate restart files for the first equilibration is as follows:
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str psffile = ../toppar/villin_crowding.psf # CHARMM protein structure file pdbfile = ../toppar/villin_crowding.pdb # PDB file reffile = ../toppar/villin_crowding.pdb rstfile = ../4_heat/run().rst [OUTPUT] rstfile = setup().rst # DCD trajectory file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NPT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat [SELECTION] group1 = an:N|an:CA|an:C|an:O [RESTRAINTS] nfunctions = 1 # number of functions function1 = POSI # restraint function type constant1 = 10.0 # force constant select_index1 = 1 # restraint group [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.6 Equilibration with restraints on the macromolecules
The control input for the first equilibration is shown below. ensemble is changed from NVT to NPT.
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str rstfile = ../5_prst_eq1/setup().rst [OUTPUT] dcdfile = run().dcd # DCD trajectory file rstfile = run().rst # restart file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] integrator = LEAP # [LEAP,VVER] nsteps = 10000 # number of MD steps timestep = 0.002 # time step (ps) eneout_period = 100 # energy output period crdout_period = 500 # coordinates output period rstout_period = 5000 # restart output period [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NPT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat temperature = 298.15 # initial temperature (K) [SELECTION] group1 = an:N|an:CA|an:C|an:O [RESTRAINTS] nfunctions = 1 # number of functions function1 = POSI # restraint function type constant1 = 10.0 # force constant select_index1 = 1 # restraint group [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.7 Making restart for equilibration 2
The control input to generate restart files for the second equilibration without the restraints is shown below.
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str psffile = ../toppar/villin_crowding.psf # CHARMM protein structure file pdbfile = ../toppar/villin_crowding.pdb # PDB file rstfile = ../6_eq1/run().rst [OUTPUT] rstfile = setup().rst # DCD trajectory file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NPT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.8 Equilibration without restraints
The control input for the second equilibration is shown below. Any restraints are removed.
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str rstfile = ../7_prst_eq2/setup().rst [OUTPUT] dcdfile = run().dcd # DCD trajectory file rstfile = run().rst # restart file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] integrator = LEAP # [LEAP,VVER] nsteps = 10000 # number of MD steps timestep = 0.002 # time step (ps) eneout_period = 100 # energy output period crdout_period = 500 # coordinates output period rstout_period = 5000 # restart output period [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NPT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat temperature = 298.15 # initial temperature (K) [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.9 Making restart for production
In this tutorial, we use NVT ensemble for production. Therefore, you should regenerate rst files because of the change of [ENSEMBLE] section. The control input to generate restart files for production is as follows:
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str psffile = ../toppar/villin_crowding.psf # CHARMM protein structure file pdbfile = ../toppar/villin_crowding.pdb # PDB file rstfile = ../8_eq2/run().rst [OUTPUT] rstfile = setup().rst # DCD trajectory file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NVT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.10 Production
During the production run, stoptr_period is newly added in the [DYNAMICS] section. This treatment is needed to remove the translational and rotational motions of the whole system.
The control input for production is as follows:
[INPUT] topfile = ../toppar/top_all36_prot.rtf # CHARMM topology file parfile = ../toppar/par_all36_prot.prm # CHARMM parameter file strfile = ../toppar/toppar_water_ions.str rstfile = ../9_prst_prod/setup().rst [OUTPUT] dcdfile = run().dcd # DCD trajectory file rstfile = run().rst # restart file [ENERGY] electrostatic = PME # [CUTOFF,PME] switchdist = 10.0 # switch distance cutoffdist = 12.0 # cutoff distance pairlistdist = 13.5 # pair-list cutoff distance pme_ngrid_x = 144 # grid size_x in [PME] pme_ngrid_y = 144 # grid size_y in [PME] pme_ngrid_z = 144 # grid size_z in [PME] [DYNAMICS] integrator = LEAP # [LEAP,VVER] nsteps = 10000 # number of MD steps timestep = 0.002 # time step (ps) eneout_period = 100 # energy output period crdout_period = 500 # coordinates output period rstout_period = 5000 # restart output period stoptr_period = 1 [CONSTRAINTS] rigid_bond = YES # constraints of all bonds involving hydrogen [ENSEMBLE] ensemble = NVT # [NVE,NVT,NPT] tpcontrol = LANGEVIN # thermostat temperature = 298.15 # initial temperature (K) [BOUNDARY] type = PBC # [PBC,NOBC] box_size_x = 150.0 # box size (x) in [PBC] box_size_y = 150.0 # box size (y) in [PBC] box_size_z = 150.0 # box size (z) in [PBC] domain_x = 2 # # of domain for x-axis domain_y = 2 # # of domain for y-axis domain_z = 2 # # of domain for z-axis
3.1.11 Convert
After the production run, you need to combine multiple trajectory files using pcrd_convert.
The control input for assembling the trajectories from production run is as follows:
[INPUT] psffile = ../toppar/villin_crowding.psf # protein structure file reffile = ../toppar/villin_crowding.pdb # PDB file [OUTPUT] pdbfile = prod.pdb # PDB file trjfile = prod.dcd # trajectory file rmsfile = prod.rms # RMSD file trrfile = prod.trr # TRROT file [TRAJECTORY] trjfile1 = ../10_prod/run().dcd # trajectory file md_step1 = 10000 # number of MD steps mdout_period1 = 500 # MD output period ana_period1 = 500 # analysis period [SELECTION] group1 = all # selection group 1 [FITTING] fitting_method = NO # method [OPTION] check_only = NO # (YES/NO) trjout_format = DCD # (PDB/DCD) trjout_type = COOR+BOX # (COOR/COOR+BOX) trjout_atom = 1 # atom group pbc_correct = NO # (NO/MOLECULE) rename_res1 = HSE HIS rename_res2 = HSD HIS
In this tutorial, no fittings of macromolecule orientations and no corrections for periodic boundary condition were applied. The assembled trajectory shows continuous diffusion from the initial position. In the [TRAJECTORY] section, you should specify the same number of md_step, and mdout_period as those in the MD simulation which generated the trajfile. ana_period can control the frequency of output in the combined trajectory. Note that pcrd_convert does not support MPI, but OpenMP.
These are the sample of command for each step using the combination of 12MPI x 8 OpenMP
# 3.1.1 Making restart for minimization $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/prst_setup prst_min.inp >prst_min_log # 3.1.2 1000 step energy minimization # with position restraint for backbone $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/spdyn min.inp > min_log # 3.1.3 Making restart for heat up $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/prst_setup prst_heat.inp >prst_heat_log # 3.1.4 5000 step heat up with position restraint under the NVT condition. $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/spdyn heat.inp > heat_log # 3.1.5 Making restart for eqauilibration with restraint $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/prst_setup prst_eq1.inp >prst_eq1_log # 3.1.6 10000 step equilibration with position restraint under the NPT condition. $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/spdyn heat.inp > heat_log # 3.1.7 Making restart for eqauilibration without restraint $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/prst_setup prst_eq2.inp >prst_eq2_log # 3.1.8 10000 step equilibration without any restraint under the NPT condition. $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/spdyn heat.inp > heat_log # 3.1.9 Making restart for production under the NVT condition $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/spdyn prst_prod.inp >prst_prod_log # 3.1.10 10000 step production without any restraint under the NVT condition. $ export OMP_NUM_THREADS=8 $ mpirun -np 8 /home/user/genesis/bin/spdyn prod.inp > prod_log # 3.1.11 combine the multiple trajectories from production run $ export OMP_NUM_THREADS=8 $ /home/user/genesis/bin/pcrd_convert pcrd_conv.inp > pcrd_conv_log
Written by Isseki Yu@RIKEN Theoretical molecular science laboratory
Oct 20, 2016.