environment/Script_run_c2_0.py at main · MMunibas/environment · GitHub

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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# Test Script to import PhysNet as energy function in CHARMM via PyCHARMM

# Basics
import os
import sys
import ctypes
import pandas
import numpy as np

# ASE
from ase import Atoms
from ase import io
import ase.units as units

# PyCHARMM
import pycharmm
import pycharmm.generate as gen
import pycharmm.ic as ic
import pycharmm.coor as coor
import pycharmm.energy as energy
import pycharmm.dynamics as dyn
import pycharmm.nbonds as nbonds
import pycharmm.minimize as minimize
import pycharmm.crystal as crystal
import pycharmm.image as image
import pycharmm.psf as psf
import pycharmm.read as read
import pycharmm.write as write
import pycharmm.settings as settings
import pycharmm.lingo as stream
import pycharmm.select as select
import pycharmm.shake as shake
import pycharmm.cons_fix as cons_fix
import pycharmm.cons_harm as cons_harm
from pycharmm.lib import charmm as libcharmm
import pycharmm.lib as lib

# Step 0: Load parameter files
#-----------------------------------------------------------

# Residue and classical force field parameter

rtf_fn = "top_all22_prot_c2.inp"
read.rtf(rtf_fn)

prm_fn = "par_all22_prot_c2.inp"
read.prm(prm_fn)

settings.set_bomb_level(-2)
settings.set_warn_level(-1)

# Step 1: Read System
#-----------------------------------------------------------

stream.charmm_script('set name c2')

stream.charmm_script('set num 0')

psf_fn = "c2.psf"
read.psf_card(psf_fn)

coor_fn = "c2.cor"
read.coor_card(coor_fn, flex=True)

# Save system pdb files
write.coor_pdb("init_c2.pdb")


# Step 2: Define PhysNet energy function
#-----------------------------------------------------------

# Prepare PhysNet input parameter
selection = pycharmm.SelectAtoms(seg_id='CH3C')

# Atomic numbers
ase_pept = io.read("init_c2.pdb", format="proteindatabank")
Z = ase_pept.get_atomic_numbers()
print(Z)

# Checkpoint files
checkpoint = "../physnet-criegee/criegee-103635_b"

# PhysNet config file
config = "../physnet-criegee/config"

# Model units are eV and Angstrom
econv = 1./(units.kcal/units.mol)
fconv = 1./(units.kcal/units.mol)

charge = 0

# Initialize PhysNet calculator
pycharmm.MLpot(
    selection,
    fq=True,
    Z=Z,
    checkpoint=checkpoint,
    config=config,
    charge=charge,
    econv=econv,
    fconv=fconv,
    v1=True)	# Model is trained by PhysNet using tensorflow 1.x

# Custom energy
energy.show()


# Step 3: Simulation - CHARMM, PhysNet
#-----------------------------------------------------------

simulation = """

energy
mini sd nstep 750
mini abnr nstep 100

write coor card name mini_@name_@num.cor

!Heating
set proc heat
open unit 11 write form name  @name_@num_@proc.res

dynamics start time 0.0002 nstep 20000 -
firstt 48.0 finalt 300.0 teminc 10.0 ihtfrq 500 -
ieqfrq 2000 ichecw 1 twindl -5.0 twindh +5.0 iasors 0 -
nprint 100 iprfrq 500 -
!atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
nbonds atom cdie cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 -
IUNREA -1 iunwrit 11 iuncrd -1 nsavc 100 ICHECW 1


!Equilibration
open unit 10 read form name    @name_@num_@proc.res
set proc eqb
open unit 11 write form name   @name_@num_@proc.res
open write unit 20 file name   @name_@num_@proc.dcd

dynamics restart VERL time 0.0001 nstep 1000000 -
firstt 300.0 finalt 300.0 -
nprint 1000 iprfrq 500 NTRFRQ 0 -
!atom cdie fshift vshift cutnb 14.0 ctofnb 12.0 -
nbonds atom cdie cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 -
iunread 10 iunwrit 11 iuncrd 20 nsavc 1000 ICHECW 1

open write unit 10 card name @name_@num_@proc.pdb
write coor unit 10 pdb
close unit 10

open write unit 10 card name @name_@num_@proc.cor
write coor unit 10 card
close unit 10


!Production
open unit 10  read form name @name_@num_@proc.res
set proc dyna
open write unit 11 file name @name_@num_@proc.dcd
open unit 13 write form name @name_@num_@proc.res

dynamics restart VERL time 0.0001 nstep 20000000 -
iunread 10 iunwrit 13 iuncrd 11 -
nprint 1000 nsavc 10 iprfrq 500 NTRFRQ 0 -
nbonds atom cdie cutnb 14.0 ctofnb 12.0 -
inbfrq -1 ihbfrq 0 IASVEL 0 ISCALE 0

open write unit 10 card name @name_@num_@proc.pdb
write coor unit 10 pdb
close unit 10

open write unit 10 card name @name_@num_@proc.cor
write coor unit 10 card
close unit 10

"""
stream.charmm_script(simulation)