QM/FQ coupled to MD sampling¶

a) Drawing the solute¶

The first thing you need to do is create the solute, in this case, a 2-Methyloxirane molecule. You could build it yourself in AMSinput for a molecule inside AMSinput, or just import it from any database. The steps are as follows:

1. Start AMSinput via SCM → New Input

2. Select the atoms-tool, add hydrogens, edit bonds, etc. You can also click into the search box () and type “methyloxirane” and select methyloxirane (ADF) from the drop-down menu. The ‘(ADF)’ string in the search results means that the molecule has already been optimized by ADF, using the BP86 XC potential with a TZP basis set and small core.

3. If you build your own structure, you can pre-optimize it by clicking on

After loading or drawing the structure, you will see it in the drawing area of the molecule editor (the dark area on the middle left side). Click somewhere in the drawing area to deselect all the atoms.

You can rotate, translate, and zoom your molecule using the mouse.

b) Creating the box¶

Now you have to put the solute inside a box, change the third dimension of the lattice vectors, and adjust them depending on the size of the solute and the kind of box you are interested in. To do that,

1. Open the builder tool by clicking: Edit → Builder

2. Make a box of 30x30x30 Angstrom by typing 30 into the diagonal elements of the Lattice vectors.

For larger solvated molecules it is recommended to use a larger box.

c) Adding the solvent¶

By using the Builder tool in AMSinput, you can easily fill the box with solvent molecules, water for this case, and know the density of the entire system after the water has been added. Again in the Builder panel:

1. Type ‘water’ in the line with ‘Fill box with’, in the field after copies of:

2. Select Water (ADF) from the search results

3. Specify 897 copies. Check the predicted density: 0.9972 g/mL

4. Click the Generate Molecules button on the bottom

Notice that the solvent molecules are generated at random positions and orientations, with the constraint that all atoms (between different molecules) are at least the specified distance (2.5 Angstrom) apart.

5. Close the Builder tool by clicking the Close button

6. Select Force Field engine by switching to the Force Field panel: →

7. Set Periodicity → Bulk. It means 3-dimensions

8. Click View → Periodic → Show Unit Cell to visualize the box

9. Rotate to your favorite view

Turn on/off periodic view: View → Periodic → Repeat Unit Cells or use the bottom to activate/deactivate the periodic view tool.

a) Configuring the duration and temperature of the equilibration trajectory¶

1. Click on next Task → Molecular Dynamics to go to the Model → MD panel

2. Configure 100000 steps with a time step of 0.5 fs (Click on a unit to change the unit, your choice will be remembered). This will result in a 0.05 ns (50 ps) long trajectory.

3. Set the sampling frequency to 100. With a total of 100000 steps, this will result in 1000 recorded samples.

4. Set the Initial temperature to 300K.

5. Set Initial velocities to Random

The time step should not be set larger than 0.5 fs because of OH and CH bonds. In an MD run with fixed OH and CH bond lengths one could use a larger time step size, like 2 fs. For some systems one may need a larger number of steps to reach equilibrium, but that is not needed in this case.

b) Defining the thermostat for equilibration¶

You have set an initial temperature of 300 K, but to maintain this temperature throughout the simulation it is also necessary to attach a thermostat to the system. In the equilibration step the Berendsen thermostat will be used.

1. Click on next to the Thermostat to go to the Model → Thermostat panel

2. Click the button to add a thermostat to the simulation

3. Select Thermostat → Berendsen, which can be used for equilibration

(note that NHC (Nosé–Hoover chain) is the preferred thermostat for production runs)

4. Set 300 K as the Temperature

5. Set 100 fs as the damping constant. Click on a unit to change the unit

It is also possible to start with a low-temperature MD to relax the initial set-up, or to define different temperatures for the components of the system, using several thermostats for different regions. To that end, you first must define two new regions by using the Regions panel, but it is out of the scope of this tutorial, so, for your solvated 2-Methyloxirane, you will use a single thermostat.

6. Click on next to MD main options

c) Defining the barostat for equilibration¶

Just like using a Thermostat to control the temperature of the system, a Barostat can be applied to keep the pressure constant by adjusting the volume.

1. Click on next to Barostat to go to the Model → Barostat panel.

2. Select Berendsen from the menu Barostat

3. Set the desired Pressure to 1.0 atm.

4. Set the Damping constant to 500 fs. Click on a unit to change the unit

5. Select XYZ from the menu Equal (such that the box stays cubic)

6. Click on next to MD main options

d) Run the MD for equilibration¶

Now you can run your setup.

1. Use the File → Run command.

2. When asked to save your input, save it with a name of your choice. Make sure you select the Tutorial directory that you made

3. The AMSjobs window comes to the front and your job starts running.

Once your job starts running, AMSjobs will show the progress of the calculation: the last few lines of the logfile:

Note that while running, the job status symbol in AMSjobs changes. If you wish to see the full logfile while the calculation is running, just click on the logfile lines in the AMSjobs window, and it will show the logfile in the AMStail window.

In the calculation first the atom types for the GAFF force field are determined, which will take some time. Next the MD steps are calculated. Let the calculation run for at least half an hour. While it is running, you can already follow its progress in AMSmovie.

4. In the AMSInput window, click SCM → Movie

The graph on the right-hand side shows the energy as a function of the trajectory step.

You can explore different MD properties along the trajectory when it is ready or even if the simulation is still running.

Use the Graph → Add Graph menu command

Click MD Properties → Temperature

Now you have two graphs. One of them is the ‘active’ graph. When you make a new graph it will always be the active graph. You can also make a graph active by clicking on it.

You can show several graphs for different properties at the same time. Temperature and Density are useful properties to evaluate the quality and convergence of your MD simulation.

Let the calculation finish. This might take around a few hours, depending on your computer resources.

Wait until AMStail shows ‘Job … has finished’ as the last line

In the dialog that pops up, Select the ‘Use MD velocities for AMS MD restart’ checkbox

click ‘Yes, new job’ to update the coordinates and make a new job

a) Defining the thermostat of the production run¶

1. Go to the Model → Thermostat panel

2. Select Thermostat → NHC (Nosé–Hoover chain), the preferred thermostat for production runs

b) Run the MD production simulation¶

Now you can run your setup.

1. Use the File → Run command.

2. When asked to save your input, save it with a name of your choice. Make sure you select the Tutorial directory that you made

3. The AMSjobs window comes to the front and your job starts running.

Let the calculation run for at least half an hour. While it is running, you can already follow its progress in AMSmovie.

Let the calculation finish. This might take around a few hours, depending on your computer resources.

a) Extracting snapshots from the MD runs¶

Once your trajectory is ready, you can open it with AMSMovie and save the geometry for any snapshot.

1. From the AMSjobs window, select the MD production job and click SCM → Movie

2. You can select your desired frame by moving the slider through the steps of the trajectory. Use the left and right arrow keys to single-step through the frames.

3. In AMSmovie, Use the File → Save geometry command

4. When asked to save your geometry, the suggested filename contains the number of the frame (ams.150.xyz) but you can save it with a name of your choice.

5. Do this for one more frames.

You can also save the xyz coordinates for the entire trajectory by using File → Export trajectory as → XYZ(.xyz) but in this tutorial, we will use just two of the snapshots.

Now you can close all windows that belong to this tutorial:

Select the SCM → Quit All command in any AMS-GUI window

AMS users should continue with the next Route # 2 to set up the QM/FQ calculations (see below).

Route # 2 (MD-experienced users and AMS users who already run the MD)

MD-experienced users can start the tutorial here if they already have at least two snapshots extracted from any trajectory.

b) Removing PBC and cutting spheres¶

This step is to be carried out if you have not removed periodic conditions yet, or if you have not cut your snapshot yet, or if you want to reduce the size of the sphere. In the AMSjobs window,

1. Start AMSjobs

2. Start AMSinput via SCM → New Input

3. In the AMSinput window, Select File → Import Coordinates… and select the file containing the data for the snapshot, it can be either a .pdb or a .xyz file extension. For the present case, it will be ‘ams.150.xyz’

4. Use a supercell if the size of the box is smaller than the size of the solute + 2 times the radius of the cutting sphere + 2 times the size of the solvent. Select Edit → Crystal → Generate Super Cell…

5. Select one of atoms of a water molecule. Select → Select Molecule. Select → Select Atoms Of Same Type. Select View → Molecule → WireFrame (or Hidden). Only the carbon atoms of the solute should be visible as balls now.

6. Click at any place to deselect all atoms

1. Select one of the carbon atoms belonging to (one of) the solute(s).

2. Select Edit → Set Origin and select Edit → Crystal → Map Atoms To (-0.5..0.5) in order to make this atom the center of the box.

3. Select ADF engine by switching to the ADF panel: →

4. Select → Select Molecule. Your entire solute molecule should appear selected.

5. Select → Select within Radius. Write a value for the distance, in Angstroms and Press OK. The radius assures the convergence of the computed data. For small solutes, we suggest cutting 12 Å sphere centered at the solutes’ geometric center, but the choice always depends on the solute size.

6. Select → Select Molecule. This step is to complete the water molecules placed at the boundaries of the sphere

7. Select → Invert Selection

8. Finally, Press backspace key on your keyboard, or go to Atoms → Delete Atoms in the GUI to remove the molecules outside the sphere.

c) Definition of the type of calculation, the different regions of the two-layer scheme, and their boundaries¶

With the sphere already cut, you can proceed with the QM/FQ settings. In this tutorial you will run UV-Vis calculations, thus, in the ‘Main’ space of the panel bar,

1. Select Task → Single Point

2. Select XC functional → Hybrid: B3LYP

3. Select Relativity → Scalar

4. Select Basis set → TZP and Core → None

5. Select Numerical Quality → Normal

Now, using the panel bar Properties

6. Click on Properties → Excitations (UV/Vis), CD command

7. For the ‘Type of excitations’ option, Select ‘SingletOnly’

Being a QM/MM calculation, we also need to define two regions: one for the solute, and one for the water. To set up these regions:

1. Click at any place to deselect all atoms

2. Select the 2-methyloxirane molecule. You can do this by selecting a single atom of this molecule and Select → Molecule

3. In the panel bar, click the Model → Regions command

4. Click the button to add a region

4. Rename Region_1 to Solute

5. Select all other atoms of the complex (by Select → Invert Selection)

6. Click the button again, and rename Region_2 to Solvent

7. Click on the right arrow at the end of the ‘Solvent’ line. Use the ‘Hidden’ command from the menu that appears. In that way, the molecules in the ‘Solvent’ region are hidden, but the ‘Solvent’ region is still visible because it is colored.

Concerning parametrization, the GUI has parameters for QM/FQ in case water is the solvent in the FQ region.

1. Use the panel bar Model → QM/FQ command

2. Check the Enable check button

3. Click the check button ‘QM part’ for the ‘Solute’ region

4. Click the check button ‘FQ part’ for the ‘Solvent’ region

The FQ parameters (\(\chi, \eta\)) will appear in the window.

This is all the setup you need. You are now ready to run the QM/FQ calculation

5. Use the File → Run command.

6. When asked to save your input, save it with a name of your choice. For example, ‘MOX_QMFQ_frame150’

7. The AMSjobs window comes to the front and your job starts running.

8. When asked for a new setup in the dialog that pops up, click ‘No’.

d) Repeat calculation for a different snapshot¶

Follow steps 4b and 4c for a different snapshot. Use a different job name, for example, ‘MOX_QMFQ_frame432’