Difference between revisions of "Template:Team:Groningen/CONTENT/MODELING/Supporting material"

Revision as of 23:53, 18 September 2015

Introduction to molecular dynamics

Molecular dynamics is used to simulate the the physical motion of atoms and molecules. It can for example be used to simulate how lipids of a cell membrane behave, how a fictitious material would perform, if a protein would denaturalize in water or how a liposome fuses with a membrane. The simulation itself is based on Newton's laws of motion, that define classical mechanics, and force fields. These laws are as follows:

If no force is applied, the object continues at it current speed.

The acceleration of the object is equal to its mass and the sum of the forces applied to the object.

When an object influences another object then the force on the second object is equal to the exerted force, but in the opposite direction.

These laws by themselves are not enough to simulate the motions of atoms and molecules. The forces that are mention in these laws have to be defined. The forces are defined by force fields in molecular dynamics simulations. Depending on the contents of the simulation some may be better than others, because one could model electrostatic interactions well while others are good for studying phase separations. Often these force fields are based on experimental data and have certain trade offs, therefore even for something as common as water different force fields exist.

There are however a few limitations to keep in mind. First of all, because these simulations are build from small elements, the simulations grows very fast. This in turn means that the simulations can take a long time even though the systems dimensions are usually on the scale of nanometers. Not only can the systems take long because they contain a lot of elements, but the effect of interest happens after a long time (the nanoseconds and microseconds range). Secondly because classical molecular dynamics is based on Newton’s laws of motion, it does not have chemical reactions. Lastly the behavior the simulations exhibit strongly depend on the input and the settings. For example the initial position might cause excessive forces and thus crashes the program.

Scripts

PGA topology creation

The following python 2.7 code can be used to generate a PGA topology file.


/usr/bin/python
namef = 'pga{}.itp'
num = 100
def outputLine(f, line):
	def p(*args):
		f.write(line.format(*args))
	return p
def printAtoms(f, num):
	f.write(';nr\ttype\tresnr\tresidu\tatom\tcgnr\tcharge\n')
	pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{0}\t{5: d}\n')
	pa(1, 'Qd', 1, 'PGATN', 'NCC', 1)
	pa(2, 'Qa', 1, 'PGATN', 'OCO', -1)
	pa(3, 'Na', 1, 'PGATN', 'CCO', 0)
	for x in range(0, num - 2):
		i = 3 * x + 4
		pa(i, 'Nd', x + 2, 'PGAM', 'NCHC', 0)
		pa(i + 1, 'Qa', x + 2, 'PGAM', 'OCO', -1)
		pa(i + 2, 'Na', x + 2, 'PGAM', 'CCO', 0)
	i = 3 * num - 2
	pa(i, 'Nd', num, 'PGATO', 'NCHC', 0)
	pa(i + 1, 'Qa', num, 'PGATO', 'OCO', -1)
	pa(i + 2, 'Qa', num, 'PGATO', 'OCO', -1)


def printBonds(f, num):
	f.write(";i\tj\tfunct\tlength\tforce.c.\n")
	pb = outputLine(f, '{0}\t{1}\t1\t{2}\t{3}\n')
	pb(1, 2, 0.23, 18000)
	pb(1, 3, 0.33, 4000)
	for x in range(0, num - 2):
		i = 3 * x + 3
		pb(i, i + 1, 0.28, 12000)
		pb(i + 1, i + 2, 0.24, 22000)
		pb(i + 1, i + 3, 0.32, 2000)
	i = 3 * num - 2
	pb(i - 1, i, 0.28, 12000)
	pb(i, i + 1, 0.24, 18000)
	pb(i, i + 2, 0.36, 4000)


def printAngles(f, num):
	f.write(';i\tj\tk\tfunct\tangle\tforce.c.\n')
	pa = outputLine(f, '{0}\t{1}\t{2}\t2\t{3}\t{4}\n')
	pa(2, 1, 3, 111.013, 100)
	pa(1, 3, 4, 141.237, 55)
	for x in range(0, num - 2):
		i = 3 * x + 3
		# next one is the same as the first one of the terminal end
		pa(i, i + 1, i + 2, 114.646, 25)
		pa(i, i + 1, i + 3, 124.643, 250)
		pa(i + 1, i + 3, i + 4, 133.794, 250)
		pa(i + 2, i + 1, i + 3, 113.393, 40)
	i = 3 * num - 3
	pa(i, i + 1, i + 2, 114.646, 25)
	pa(i, i + 1, i + 3, 127.488, 100)
	pa(i + 2, i + 1, i + 3, 116.572, 45)
def main():
	with open(namef.format(num), 'w') as f:
		f.write('[ moleculetype ]\n;molname\tnrexcl\nPGA{}\t\t1\n'.format(num))
		f.write('\n[ atoms ]\n')
		printAtoms(f, num)
		f.write('\n[ bonds ]\n')
		printBonds(f, num)
		f.write('\n[ angles ]\n')
		printAngles(f, num)
if __name__ == "__main__":

   main()