monte carlo simulation of protein structure and grid

Question

I am working Monte Carlo simulation script over protein structure. I have never done before Monte Carlo scripting. I will extent this program at large scale. According to protein xyz coordinates I have to define the box size. This box will be divided into a grid of size 0.5 A. Based on distance and angle criteria I have to assign the point based on Boltzmann probability distribution.

My program should be move in each direction by taking grid of 0.5 A and generate the random point and check the condition of distance and angle. If satisfy the condition put point there otherwise discard that point based on Boltzmann probability distribution.

Here is my code for generation of random points

from __future__ import division    
import math as mean    
from numpy import *   
import numpy as np   
from string import *    
from random import *    

def euDist(cd1, cd2):# calculate distance
    d2 = ((cd1[0]-cd2[0])**2 + (cd1[1]-cd2[1])**2 + (cd1[2]-cd2[2])**2)
    d1 = d2 ** 0.5
    return round(d1, 2)

def euvector(c2,c1):# generate vector
    x_vec = (c2[0] - c1[0])
    y_vec = (c2[1] - c1[1])
    z_vec = (c2[2] - c1[2])
    return (x_vec, y_vec, z_vec)


 for arang in range(1000):  # generate random point
        arang = arang + 1
        x,y,z = uacoord
        #print x,y,z

        x1,y1,z1 = (uniform(x-3.5,x+3.5), uniform(y-3.5,y+3.5), uniform(z-3.5,z+5))
        pacord = [x1,y1,z1]                 # random point coordinates
        print pacord

I am completely struck to generate the box size from the xyz coordinates of protein structure and how to define the grid of size 0.5. How to check every point in the box.
Any help will be appreciable.

Answer 1

I love your topic, question, and approach. I don't know how long it'll stick around here.

A 3D mesh in a rectangular space is common in finite element analysis and all other techniques for analyzing physics problems. Have a look at how meshes are generated.

There are usually two parts: a mesh of (x,y,z) points in space and boxes that connect them. A simple 2D mesh with two elements would look like this:

D               E               F
o (1, 0) ------ o (1, 1) ------ o (1, 2)
+               +               +
+               +               +
+               +               +
o (0, 0) -------+ (0, 1) -------+ (0, 2)
A               B               C

There are six points in this mesh:

A (0, 0)
B (0, 1)
C (0, 2)
D (1, 0)
E (1, 1)
F (1, 2)

and two boxes:

1 - (A, B, E, D)
2 - (B, C, F, E)

So one possible approach would be to iterate over all the boxes, check position at the centroid, and adjust accordingly.

I'd externalize the mesh definition from your code and read it in from a file. That way you can handle different meshes with the same code.

If I understand you correctly, the mesh will remain fixed in space. You're trying to randomize the motion of the protein itself. It's like a fluids problem, an Eulerian approach: you're tracking the motion of material against a fixed control volume. Except it's a protein instead of a fluid.

So you'll also have a separate definition of the initial condition for the protein in space. You plan to increment it in time to see how the protein folds according to your rules.

Am I close?

Answer 2

My code was written for a protein folding application, but the overall idea is the same. It starts with a certain temperature and decreases the temperature step wise, the protein (or in your case the 'point) is moved randomly, if the energy score of the new position/structure is lower than the previous one it is accepted, if not the pose will be evaluated according to the Metropolis distribution. You have to define your scorefunction(), a function to define a random start position and a mover which moves your point from its original position. The code below is modified from my original script, just to give you a rough idea.

kT_lower=0.1
kT_upper=100
ktemp=kT_upper

max_iterations=15000

i=-1

#create random start point
pose=create_pose()

#evaluate start point
starting_score=scorefunction(pose)

while i<max_iterations:
    i=i+1
    new_pose=random_move(pose)
    if scorefunction(new_pose)<scorefunction(pose):
        pose=new_pose
    else:
        #linear decrease of kT
        ktemp=kT_upper-i*(kT_upper-kT_lower)/max_iterations
        #exponentatial decrease of kT
        #ktemp=math.exp(float(i)/float(max_iterations)*float(-5))*kT_upper+kT_lower

        try:
            p=math.exp(DeltaE/ktemp)
        except OverflowError:
            p=-1

        if random.random()<p:
            pose=new_pose
            print str(i)+'; accept new pose, metropolis'
        else:
            print str(i)+'; reject new pose!'