Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • nick
  • pinyili2
  • dev
  • 1.0a
  • 0.5a
  • 0.4
  • 0.3
  • 0.2
  • 0.1
  • 0.1.dev0
11 results
Blame Permalink
simulate.py 5.61 KiB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131 






import os
import tempfile
from .coords import readArbdCoords, readAvgArbdCoords
import shutil
from . import get_resource_path

arbd="/home/cmaffeo2/development/cuda/arbd.dbg/src/arbd" # reduced the mem footprint cause vmd
namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multicore-CUDA/namd2"

def multiresolution_simulation( model, output_name, 
                                job_id=None,
                                gpu = 0,
                                directory=None,
                                coarse_steps = 5e7,
                                fine_steps = 5e7,
                                coarse_output_period = 1e5,
                                fine_output_period = 1e5,
                                remove_long_bonds=False,
                                backbone_scale=1.0,
                            ):

    ## Round steps up to nearest multiple of output_period, plus 1
    coarse_steps = ((coarse_steps+coarse_output_period-1)//coarse_output_period)*coarse_output_period+1
    fine_steps = ((fine_steps+fine_output_period-1)//fine_output_period)
    if fine_steps == 1: fine_steps += 1
    fine_steps = fine_steps*fine_output_period+1 

    ret = None
    d_orig = os.getcwd()
    try:
        if directory is None:
            tmp = job_id if job_id is not None else output_name
            directory = tempfile.mkdtemp(prefix='origami-%s-' % tmp, dir='/dev/shm/')
        elif not os.path.exists(directory):
            os.makedirs(directory)
        os.chdir(directory)

        output_directory = "output"
        
        # TODO: add the output directory in to the readArbdCoords functions or make it an attribute of the model object 

        """ Coarse simulation """
        simargs = dict(timestep=200e-6, output_period=coarse_output_period, gpu=gpu)

        model._clear_beads()
        model._generate_bead_model( 5, 5, local_twist=False )

        if remove_long_bonds:
            output_prefix = "%s-0" % output_name
            full_output_prefix = "%s/%s" % (output_directory,output_prefix)
            ## TODO Remove long bonds
            model.simulate( output_name = output_prefix, numSteps=0.05*coarse_steps, **simargs )
            try:
                coordinates = readArbdCoords('%s.restart' % full_output_prefix)
            except:
                coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
            output_prefix = "%s-1" % output_name
            model.update_splines(coordinates)
            model.simulate( output_prefix = output_prefix, num_steps=0.95*coarse_steps, **simargs )
        else:
            output_prefix = "%s-1" % output_name
            full_output_prefix = "%s/%s" % (output_directory,output_prefix)
            model.simulate( output_name = output_prefix, num_steps=coarse_steps, **simargs )
        try:
            coordinates = readArbdCoords('%s.restart' % full_output_prefix)
        except:
            coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)


        """ Fine simulation """ 
        output_prefix = "%s-2" % output_name
        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
        simargs['timestep'] = 50e-6
        simargs['output_period'] = fine_output_period
        model.update_splines(coordinates)
        model._clear_beads()
        model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
        model.simulate( output_name = output_prefix, num_steps=fine_steps, **simargs )
        if fine_steps <= 2*fine_output_period+1:
            try:
                coordinates = readArbdCoords('%s.restart' % full_output_prefix)
            except:
                coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
        else:
            try:
                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
            except:
                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)


        """ Freeze twist """
        output_prefix = "%s-3" % output_name
        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
        model.update_splines(coordinates)
        model._clear_beads()
        model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=False )
        model.simulate( output_name = output_prefix, num_steps=fine_steps, **simargs )
        if fine_steps <= 2*fine_output_period+1:
            try:
                coordinates = readArbdCoords('%s.restart' % full_output_prefix)
            except:
                coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
        else:
            try:
                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
            except:
                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)


        """ Atomic simulation """
        output_prefix = "%s-4" % output_name
        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
        model.update_splines(coordinates)
        model._clear_beads()
        model._generate_atomic_model(scale=backbone_scale)
        try:
            shutil.copytree( get_resource_path("charmm36.nbfix"), "charmm36.nbfix" )
        except FileExistsError:
            ...
        model.write_atomic_ENM( output_prefix )
        model.atomic_simulate( output_name = output_prefix )

        ret = directory
        
    except:
        raise
    finally:
        os.chdir(d_orig)
        
    return ret