一、 pymol处理3LFT

二、3LFT_apo_GaMD  

清理pdb文件（clean up Protein PDB files）

修改部分残基名

2.1  pdb清理

2.2 质子化

conda activate pdb2pqr
pdb2pqr30[hX8]  3LFT_apo_clean2.pdb 3LFT_apo_clean2.pqr --ff AMBER --ffout AMBER --with-ph 7.4 --pdb-output 3LFT_apo_pqr.pdb
得到3LFT_apo_pqr.pdb
 
检查
ls -lh 3LFT_apo_clean2.pqr 3LFT_apo_pqr.pdb
grep '^ATOM' 3LFT_apo_pqr.pdb | awk '{print $4}' | sort | uniq | grep -E 'HIS|HID|HIE|HIP|CYS|CYX|LYN|ASH|GLH'
grep '^SSBOND' 3LFT_apo_pqr.pdb
grep '^ATOM' 3LFT_apo_pqr.pdb | head
grep '^ATOM' 3LFT_apo_pqr.pdb | tail
 
·  只有 HIE，说明组氨酸已经被指定为 ε 位质子化形式[zz9]  
·  没有 SSBOND 记录[zz10]  
·  没看到 CYX / ASH / GLH / HIP 这些特殊残基 
·  有正常的 OXT 末端
 
 
检查：pymol显示组氨酸
select his_all, resn HIS+HID+HIE+HIP
show sticks, his_all
color yellow, his_all
  

2.3 Tleap

cat > leap_check.in <<'EOF'
source oldff/leaprc.ff14SB
source leaprc.water.tip3p
loadAmberParams frcmod.ionsjc_tip3p
pro = loadpdb 3LFT_apo_pqr.pdb
check pro
charge pro
desc pro
quit
EOF
Warning
 
/home/chpeng/software/amber18/bin/teLeap: Warning!
The unperturbed charge of the unit (-6.000000) is not zero.
# ·  蛋白总电荷 = -6 后面要加离子中和
/home/chpeng/software/amber18/bin/teLeap: Warning!
Close contact of 1.455130 angstroms between .R<ASP 93>.A<OD2 10> and .R<THR 149>.A<HG1 12>
 
/home/chpeng/software/amber18/bin/teLeap: Warning!
Close contact of 1.159596 angstroms between .R<SER 212>.A<HG 9> and .R<LYS 214>.A<HG2 9>
 
/home/chpeng/software/amber18/bin/teLeap: Warning!
Close contact of 1.470285 angstroms between .R<SER 231>.A<HG 9> and .R<VAL 232>.A<H 2>
 

2.4 生成拓扑

cat > leap_volume.in <<'EOF'
source oldff/leaprc.ff14SB
source leaprc.water.tip3p
loadAmberParams frcmod.ionsjc_tip3p
 
pro = loadpdb 3LFT_apo_pqr.pdb
solvatebox pro TIP3PBOX 12.0
charge pro
quit
EOF
tleap -f leap_volume.in > leap_volume.log
grep -n "Volume:" leap_volume.log 52: Volume: 477561.964 A^3 (pdb2pqr) [chpeng@localhost 3LFT_apo_gamd]$ 
grep -n "Added" leap_volume.log | tail 54: Added 11364 residues. (pdb2pqr) [chpeng@localhost 3LFT_apo_gamd]$ 
grep -n "Total unperturbed charge" leap_volume.log 55:Total unperturbed charge: -6.000000
 
NaCl 盐对数 = 43 对[LW11] 
addionsrand pro Na+ 49 Cl- 43
溶剂化加盒子
cat > leap_solvate.in <<'EOF'
source oldff/leaprc.ff14SB
source leaprc.water.tip3p
loadAmberParams frcmod.ionsjc_tip3p
 
pro = loadpdb 3LFT_apo_pqr.pdb
 
solvatebox pro TIP3PBOX 12.0[zz12] [LW13] 
charge pro
 
addionsrand pro Na+ 6 
charge pro
 
saveamberparm pro pro-sol.prmtop pro-sol.inpcrd[zz14] [LW15] 
savepdb pro pro.pdb
 
quit
EOF
得到pro-sol.prmtop  pro-sol.inpcrd  pro-sol.pdb

2.5 pre-equrlibrium

01-min.in[hX16] [hX17] 
minimization simulation
&cntrl
 imin=1,       ! Minimize the initial structure
 maxcyc=1000, ! Maximum number of cycles for minimization
 ncyc=500,    ! Switch from steepest descent to conjugate gradient minimization after ncyc cycles
 ntb=1,        ! Constant volume
 ntp=0,        ! No pressure scaling
 ntf=1,        ! Complete force evaluation
 ntc=1,        ! No SHAKE
 ntpr=5000,      ! Print to mdout every ntpr steps
 ntwr=200,    ! Write a restart file every ntwr steps
 cut=9.0,      ! nonbonded cutoff, in angstroms
 
 ! Wrap coordinates when printing them to the same unit cell
 iwrap=1,
 ntwprt = 4383,#默认输出全原子[LW18] 
 ntr=1, restraintmask='(@1-4383)',（限制为全原子[zz19] [LW20] ，与cmd教程中的:1-439'    (这里的1-439修改成你的体系复合物实际的残基数目编号)  不一样
 restraint_wt=1.0
 /[浩夏21] 
02-min.in
minimization simulation without restrain
&cntrl
 imin=1,
 maxcyc=100000,
 ncyc=50000,
 ntb=1,
 ntp=0,
 ntf=1,
 ntc=1,
 ntpr=5000,
 ntwr=5000,
 cut=9.0,
 
 iwrap=1,
 ntwprt=,
/
03-nvt.in
heat the system to 310 for 200ps and then equilbriation for 200ps in nvt
&cntrl
 imin=0,irest=0,ntx=1,
 nstlim=2500,dt=0.001,
 ntc=2,ntf=2,
 ntb=1,
 cut=9.0,      ! nonbonded cutoff, in angstroms
 ntpr=500, ntwx=500,
 ntt=3, gamma_ln=1.0,
 ntr=1, restraintmask='!:WAT,Na+,K+,Cl-',（模板这里和之前的不一致[zz22] ）
 restraint_wt=0.5,
 nmropt=1,
 ! Wrap coordinates when printing them to the same unit cell
 iwrap=1,
 ntwprt = ,
 /
 &wt TYPE='TEMP0', istep1=0, istep2=2500,
   value1=0.1, value2=100.0, /(此处100)[浩夏23] [LW24] [浩伦25] 
 &wt TYPE='END'/
04-nvt.in
density the system for another 200ps
&cntrl
 imin=0,        ! No minimization
 irest=1,       ! This IS a new MD simulation
 ntx=5,         ! read coordinates only
 ! Temperature control
  ntt=3,         ! Langevin dynamics
  gamma_ln=1.0,  ! Friction coefficient (ps^-1)
  ig=-1,         ! random seed
  ! Potential energy control
  cut=9.0,       ! nonbonded cutoff, in Angstroms
 
  ! MD settings
  nstlim=50000, ! simulation length
  dt=0.002,      ! time step (ps)
  ! SHAKE
  ntc=2,         ! Constrain bonds containing hydrogen
  ntf=2,         ! Do not calculate forces of bonds containing hydrogen
  ! Control how often information is printed
  ntpr=1000,     ! Print energies every 1000 steps
  ntwx=1000,     ! Print coordinates every 1000 steps to the trajectory
  ntwr=5000,    ! Print a restart file every 10K steps (can be less frequent)
  ntxo=2,        ! Write NetCDF format
  ioutfm=1,      ! Write NetCDF format (always do this!)
  ! Wrap coordinates when printing them to the same unit cell
  iwrap=1,
  ntwprt = 4383,
  ! Constant pressure control. Note that ntp=3 requires barostat=1
  barostat=2,    ! Berendsen... change to 2 for MC barostat
  ntp=1,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/ [浩伦26] surften
（ntp=3,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/ surften
  csurften=3, gamma_ten=0.0, ninterface=2,（此处原模版为膜体系，进行更改）
  pres0=1.0,     ! Target external pressure, in bar
  taup=0.5,      ! Berendsen coupling constant (ps)
  ntr=1, restraintmask='!:WAT,Na+,K+,Cl-',
  restraint_wt=0.5,
  /
 &wt type='TEMP0', istep1=0, istep2=50000,
                    value1=100.0, value2=310.0 /[浩伦27] 
 &wt TYPE='END' /
 
05-bb
density the system for  100ps
&cntrl
 imin=0,        ! No minimization
 irest=1,       ! This IS a new MD simulation
 ntx=5,         ! read coordinates only
 ! Temperature control
  ntt=3,         ! Langevin dynamics
  gamma_ln=1.0,  ! Friction coefficient (ps^-1)
  tempi=310.0,   ! Initial temperature
  temp0=310.0,   ! Target temperature
  ig=-1,         ! random seed
  ! Potential energy control
  cut=9.0,       ! nonbonded cutoff, in Angstroms
 
  ! MD settings
  nstlim=500000, ! simulation length
  dt=0.002,      ! time step (ps)
  ! SHAKE
  ntc=2,         ! Constrain bonds containing hydrogen
  ntf=2,         ! Do not calculate forces of bonds containing hydrogen
  ! Control how often information is printed
  ntpr=1000,     ! Print energies every 1000 steps
  ntwx=1000,     ! Print coordinates every 1000 steps to the trajectory
  ntwr=5000,    ! Print a restart file every 10K steps (can be less frequent)
  ntxo=2,        ! Write NetCDF format
  ioutfm=1,      ! Write NetCDF format (always do this!)
  ! Wrap coordinates when printing them to the same unit cell
  iwrap=1,
  ntwprt = 4383,
  ! Constant pressure control. Note that ntp=3 requires barostat=1
  barostat=2,    ! Berendsen... change to 2 for MC barostat
  ntp=1,
 
  （ntp=3,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/surften
csurften=3, gamma_ten=0.0, ninterface=2,）删除
  pres0=1.0,     ! Target external pressure, in ba[浩夏28] r
  taup=0.5,      ! Berendsen coupling constant (ps)
ntr=1, restraintmask='@CA,C,N',
  ntr=1, restraintmask=':MOL | @CA,C,N',空蛋白MOL去掉[浩伦29] 
  restraint_wt=0.5,
  /
06-ca
density the system for 100ps
&cntrl
 imin=0,        ! No minimization
 irest=1,       ! This IS a new MD simulation
 ntx=5,         ! read coordinates only
 ! Temperature control
  ntt=3,         ! Langevin dynamics
  gamma_ln=1.0,  ! Friction coefficient (ps^-1)
  tempi=310.0,   ! Initial temperature
  temp0=310.0,   ! Target temperature
  ig=-1,         ! random seed
  ! Potential energy control
  cut=9.0,       ! nonbonded cutoff, in Angstroms
 
  ! MD settings
  nstlim=500000, ! simulation length
  dt=0.002,      ! time step (ps)
  ! SHAKE
  ntc=2,         ! Constrain bonds containing hydrogen
  ntf=2,         ! Do not calculate forces of bonds containing hydrogen
  ! Control how often information is printed
  ntpr=1000,     ! Print energies every 1000 steps
  ntwx=1000,     ! Print coordinates every 1000 steps to the trajectory
  ntwr=5000,    ! Print a restart file every 10K steps (can be less frequent)
  ntxo=2,        ! Write NetCDF format
  ioutfm=1,      ! Write NetCDF format (always do this!)
  ! Wrap coordinates when printing them to the same unit cell
  iwrap=1,
  ntwprt = 4383,
  ! Constant pressure control. Note that ntp=3 requires barostat=1
  barostat=2,    ! Berendsen... change to 2 for MC barostat
  ntp=1,
（ntp=3,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/surften
csurften=3, gamma_ten=0.0, ninterface=2,）删除
[浩夏30] 
  pres0=1.0,     ! Target external pressure, in bar
  taup=0.5,      ! Berendsen coupling constant (ps)
  ntr=1, restraintmask='@CA',
ntr=1, restraintmask=':MOL | @CA',空蛋白MOL去掉
  restraint_wt=0.5,
/
Md-10ns.in 进行10ns常规md
density the system for another 10 ns
&cntrl
 imin=0,        ! No minimization
 irest=0,       ! This IS a new MD simulation
 ntx=1,         ! read coordinates only
 ! Temperature control
 ntt=3,         ! Langevin dynamics
 gamma_ln=1.0,  ! Friction coefficient (ps^-1)
 tempi=310.0,   ! Initial temperature
 temp0=310.0,   ! Target temperature
 ig=-1,         ! random seed
! Potential energy control
 cut=9.0,      ! nonbonded cutoff, in angstroms
 
! MD settings
 nstlim=5000000, ! simulation length
 dt=0.002,      ! time step (ps)
 ! SHAKE
  ntc=2,         ! Constrain bonds containing hydrogen
  ntf=2,         ! Do not calculate forces of bonds containing hydrogen
  ! Control how often information is printed
  ntpr=1000,     ! Print energies every 1000 steps
  ntwx=1000,     ! Print coordinates every 1000 steps to the trajectory
  ntwr=5000,    ! Print a restart file every 10K steps (can be less frequent)
  ntxo=2,        ! Write NetCDF format
  ioutfm=1,      ! Write NetCDF format (always do this!)
  ! Wrap coordinates when printing them to the same unit cell
  iwrap=1,
  ntwprt = 4383,
  ! Constant pressure control. Note that ntp=3 requires barostat=1
  barostat=2,
  ntp=1,
ntp=3,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/surften
csurften=3, gamma_ten=0.0, ninterface=2,）删除
 
  pres0=1.0,     ! Target external pressure, in bar
  taup=0.5,      ! Berendsen coupling constant (ps)
/
 
md.in
GaMD simulation
 &cntrl
    imin=0,        ! No minimization
    irest=0,       ! This IS a new MD simulation
    ntx=1,         ! read coordinates only
 
    ! Temperature control
    ntt=3,         ! Langevin dynamics
    gamma_ln=1.0,  ! Friction coefficient (ps^-1)
    tempi=310.0,   ! Initial temperature
    temp0=310.0,   ! Target temperature
    ig=-1,         ! random seed
 
    ! Potential energy control
    cut=9.0,      ! nonbonded cutoff, in angstroms
    !fswitch=10.0,  ! Force-based switching
 
    ! MD settings
    nstlim=500000, ! simulation length
    dt=0.002,      ! time step (ps)
 
    ! SHAKE
    ntc=2,         ! Constrain bonds containing hydrogen
    ntf=2,         ! Do not calculate forces of bonds containing hydrogen
 
    ! Control how often information is printed
    ntpr= 1000,     ! Print energies every 1000 steps
    ntwx= 1000,     ! Print coordinates every 1000 steps to the trajectory
    ntwr=5000,    ! Print a restart file every 10K steps (can be less frequent)
!   ntwv=-1,       ! Uncomment to also print velocities to trajectory
!   ntwf=-1,       ! Uncomment to also print forces to trajectory
    ntxo=2,        ! Write NetCDF format
    ioutfm=1,      ! Write NetCDF format (always do this!)
 
    ! Wrap coordinates when printing them to the same unit cell
    iwrap=1,
    ntwprt = ,（与原模版31876近似。下面参数设置相同）
 
    ! Constant pressure control. Note that ntp=3 requires barostat=1
    barostat=2,    ! Berendsen... change to 2 for MC barostat
    ntp=3,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/ surften
    pres0=1.0,     ! Target external pressure, in bar
    taup=0.5,      ! Berendsen coupling constant (ps)
 
    ! Constant surface tension (needed for semi-isotropic scaling). Uncomment
    ! for this feature. csurften must be nonzero if ntp=3 above
    csurften=3,    ! Interfaces in 1=yz plane, 2=xz plane, 3=xy plane
    gamma_ten=0.0, ! Surface tension (dyne/cm). 0 gives pure semi-iso scaling
    ninterface=2,  ! Number of interfaces (2 for bilayer)
 
    ! GaMD parameters
    igamd = 3, iE = 1, irest_gamd = 0,
    ,
ntcmd = 176000, nteb = 176000, ntave = 17600,
ntcmdprep = 8800, ntebprep = 8800,[浩夏31] 
（原子数38490）
ntcmd = 11200000, nteb = 11200000, ntave = 1120000,（1106264）
    ntcmdprep = 560000, ntebprep = 560000,（553132）[浩伦32] 
    sigma0P = 12.0, sigma0D = 12.0,
 /
 
 
 
 
 
gamd-restart.in
GaMD simulation
 &cntrl
    imin=0,        ! No minimization
    irest=1,       ! This IS a new MD simulation
    ntx=5,         ! read coordinates only
 
    ! Temperature control
    ntt=3,         ! Langevin dynamics
    gamma_ln=1.0,  ! Friction coefficient (ps^-1)
    tempi=310.0,   ! Initial temperature
    temp0=310.0,   ! Target temperature
    ig=-1,         ! random seed
 
    ! Potential energy control
    cut=9.0,      ! nonbonded cutoff, in angstroms
    !fswitch=10.0,  ! Force-based switching
 
    ! MD settings
    nstlim=500000, ! simulation length
    dt=0.002,      ! time step (ps)
 
    ! SHAKE
    ntc=2,         ! Constrain bonds containing hydrogen
    ntf=2,         ! Do not calculate forces of bonds containing hydrogen
 
    ! Control how often information is printed
    ntpr=,     ! Print energies every 1000 steps
    ntwx=,     ! Print coordinates every 1000 steps to the trajectory
    ntwr=5000,    ! Print a restart file every 10K steps (can be less frequent)
!   ntwv=-1,       ! Uncomment to also print velocities to trajectory
!   ntwf=-1,       ! Uncomment to also print forces to trajectory
    ntxo=2,        ! Write NetCDF format
    ioutfm=1,      ! Write NetCDF format (always do this!)
 
    ! Wrap coordinates when printing them to the same unit cell
    iwrap=1,
    ntwprt = 3,
 
    ! Constant pressure control. Note that ntp=3 requires barostat=1
    barostat=2,    ! Berendsen... change to 2 for MC barostat
    ntp=3,         ! 1=isotropic, 2=anisotropic, 3=semi-isotropic w/ surften
    pres0=1.0,     ! Target external pressure, in bar
    taup=0.5,      ! Berendsen coupling constant (ps)
 
    ! Constant surface tension (needed for semi-isotropic scaling). Uncomment
    ! for this feature. csurften must be nonzero if ntp=3 above
    csurften=3,    ! Interfaces in 1=yz plane, 2=xz plane, 3=xy plane
    gamma_ten=0.0, ! Surface tension (dyne/cm). 0 gives pure semi-iso scaling
    ninterface=2,  ! Number of interfaces (2 for bilayer)
 
    ! GaMD parameters
    igamd = 3, iE = 1, irest_gamd = 1,   
    ntcmd = 0, nteb = 0, ntave = 17600(154000, , ntave = 1120000,（正式运行改成）
    ntcmdprep = 0, ntebprep = 0,
    sigma0P = 6.0, sigma0D = 6.0,[浩夏33] 
 /
run-Gamd.sh
 
 
#!/bin/sh
export CUDA_VISIBLE_DEVICES=1
 
export CUDA_HOME=/usr/local/cuda-9.1
export PATH=/usr/local/cuda-9.1/bin:$PATH
export LD_LIBRARY_PATH=/usr/local/cuda-9.1/lib64:/home/chpeng/software/amber18/lib:/home/chpeng/intel/compilers_and_libraries_2017.8.262/linux/compiler/lib/intel64:/home/chpeng/intel/compilers_and_libraries_2017.8.262/linux/compiler/lib/intel64_lin:/home/chpeng/intel/compilers_and_libraries_2017.8.262/linux/mkl/lib/intel64_lin:/home/chpeng/intel/compilers_and_libraries_2017.8.262/linux/mpi/intel64/lib
pmemd -O -i 01-min.in -p complex.prmtop -c complex.inpcrd -o 01-min.out -r 01-min.rst -x 01-min.nc -ref complex.inpcrd
pmemd.cuda -O -i 02-min.in -p complex.prmtop -c 01-min.rst -o 02-min.out -r 02-min.rst -x 02-min.nc -ref 01-min.rst
pmemd.cuda -O -i 03-nvt.in -p complex.prmtop -c 02-min.rst -o 03-nvt.out -r 03-nvt.rst -x 03-nvt.nc -ref 02-min.rst
pmemd.cuda -O -i 04-npt.in -p complex.prmtop -c 03-nvt.rst -o 04-npt.out -r 04-npt.rst -x 04-npt.nc -ref 03-nvt.rst
pmemd.cuda -O -i 05-bb.in -p complex.prmtop -c 04-npt.rst -o 05-bb.out -r 05-bb.rst -x 05-bb.nc -ref 04-npt.rst
pmemd.cuda -O -i 06-ca.in -p complex.prmtop -c 05-bb.rst -o 06-ca.out -r 06-ca.rst -x 06-ca.nc -ref 05-bb.rst
pmemd.cuda -O -i md-1ns.in -p complex.prmtop -c 06-ca.rst -o md-1ns.out -r md-1ns.rst -x md-1ns.nc -ref 06-ca.rst
pmemd.cuda -O -i md.in -p complex.prmtop -c md-10ns.rst -o md-200ns.out -r md-200ns.rst -x md-200ns.nc -ref md-10ns.rst
##gamd equilibration
pmemd.cuda -O -i md.in -p complex.prmtop -c md-10ns.rst -o gamd-1.out -x gamd-1.nc -r gamd-1.rst -gamd gamd-1.log -ref md-10ns.rst[浩夏34] 
##gamd production
pmemd.cuda -O -i gamd-restart.in -p complex.prmtop -c gamd-1.rst -o gamd-2.out -x gamd-2.nc -r gamd-2.rst -gamd gamd-2.log -ref gamd-1.rst
 
启动
chmod +x run-GaMD.sh
 
nohup ./run-GaMD.sh > run-GaMD.log 2>&1 &

三、3LFT_holo_GaMD    (85.9) /home/databank1/xhl/3LFT_holo_gamd/

3.1 分开蛋白配体

grep '^HETATM' 3LFT_holo.pdb | grep ' TRP A   1 ' > ligand_trp_raw.pdb
grep -v ' TRP A   1 ' 3LFT_holo.pdb > protein_holo_raw.pdb
 

3.2 TRP 配体准备[hX36] 

antechamber -i lig.sdf -fi sdf -o lig.gjf -fo gcrt -at gaff2 -gn"%nproc=4" -gm "%mem=4GB" -gk "#B3LYP/6-31G* em=gd3bj pop=MKiop(6/33=2,6/42=6) opt" -rn MOL -nc 012. nohup g16 lig.gjf
 
Maestro中ligandprep 配体预处[浩伦37] 理ph设定为7.4，输出为两个构象，区别不大。默认选择第一个。输出为lig.sdf
antechamber -i lig.sdf -fi sdf -o lig.gjf -fo gcrt -at gaff2 -gn "%nproc=4" -gm "%mem=4GB" -gk "#B3LYP/6-31G* em=gd3bj pop=MKiop(6/33=2,6/42=6) opt" -rn MOL -nc 0
 
nohup g16[浩伦38]  lig.gjf &
得到lig.log
antechamber -i lig.log -fi gout -o lig.prep -fo prepi -c resp
parmchk2 -i lig.prep -f prepi -o lig.frcmod -a y
得到lig.prep、PDB.PDB、
 

3.3 复合物准备

tleap
source leaprc.protein.ff14SB
source leaprc.gaff2
source leaprc.water.tip3p
loadamberparams lig.frcmod
loadamberprep lig.prep
mol = loadpdb lig.pdb
pro = loadpdb 3LFT_holo_pqr.pdb
com = combine {pro mol}
 
saveamberparm pro pro.prmtop pro.inpcrd
saveamberparm mol lig.prmtop lig.inpcrd
saveamberparm com native.prmtop native.inpcrd
savepdb com com-dry.pdb
 
 
检查蛋白com-dry.pdb与原来的3LFT_holo_pqr叠合
#warning不用管
charge com
 
solvatebox com TIP3PBOX 12
solvatebox com OPCBOX 12（前面加载的是 TIP3P 水模型参数）
addionsrand com Na+ 6
charge com   #为0
saveamberparm com complex.prmtop complex.inpcrd
savepdb com com.pdb
 
 
quit
得到com.pdb[浩伦39]   complex.inpcrd  complex.prmtop
查看除去水和离子的溶质数
 
 
模拟将ntwprt = 4410[LW40] 
  
Restraintmask=':MOL加上mol
(
上述一致。

四、3LFT_holo_cMD    (85.9) /home/databank1/xhl/3LFT_holo_cmd/

五、3LFT_apo_cMD    (85.9) /home/databank1/xhl/3LFT_apo_cmd/

cMD-amber.sh[浩夏41] 
 run-min.sh
 run-cMD.sh
 
85.9的vmd不显示
75.5的vmd显示
 
 
ERROR: Calculation halted.  Periodic box dimensions have changed too much from their initial values.
  Your system density has likely changed by a large amount, probably from
  starting the simulation from a structure a long way from equilibrium.
 
  [Although this error can also occur if the simulation has blown up for some reason]
 
  The GPU code does not automatically reorganize grid cells and thus you
  will need to restart the calculation from the previous restart file.
  This will generate new grid cells and allow the calculation to continue.
  It may be necessary to repeat this restarting multiple times if your system
  is a long way from an equilibrated density.
 
  Alternatively you can run with the CPU code until the density has converged
  and then switch back to the GPU code.[浩夏42] 
·  1-min.in：最小化 1000 步，不是 ns，不能直接折算真实时间 
·  02-min.in：最小化 100000 步，不是 ns 
·  03-nvt.in：nstlim=2500, dt=0.001，实际只有 2.5 ps 
·  04-npt.in：50000 × 0.002 ps = 100 ps 
·  05-bb.in：500000 × 0.002 ps = 1000 ps = 1 ns 
·  06-ca.in：500000 × 0.002 ps = 1000 ps = 1 ns 
·  md-10ns.in：5000000 × 0.002 ps = 10000 ps = 10 ns 
·  md.in：25000000 × 0.002 ps = 50000 ps = 50 ns 
·  md-continue.in：也是 50 ns

六，分析

cpptraj pro-sol.prmtop <<EOF
parmstrip :WAT,Na+,Cl-
parmwrite out pro.prmtop
run
EOF

rmsd .in

parm pro.prmtop
 
# 读取多条轨迹
trajin cmd-1.nc
trajin cmd-2.nc
trajin cmd-3.nc
trajin cmd-4.nc
 
# 处理 PBC 和居中蛋白
autoimage
center :1-4383 mass origin   # 以蛋白质残基 1-4383 为中心
 
# Backbone RMSD（Cα、C、N）以第一帧为参考结构
rms first :1-4383@CA,C,N out apo_backbone_rmsd.dat
 
run
quit

rmsf.in [hX43] 

parm pro.prmtop
 
# 读取多条轨迹
trajin cmd-1.nc
trajin cmd-2.nc
trajin cmd-3.nc
trajin cmd-4.nc
 
# PBC 处理和蛋白质居中
autoimage
center :1-4383 mass origin   # 以蛋白质残基 1-4383 为中心
 
# 对齐结构（RMS Fit），以第一帧的 CA 为参考
rms first :1-4383@CA
 
# 计算每个残基 Backbone 的 RMSF
atomicfluct :1-4383@CA,C,N byres out apo_backbone_rmsf.dat
 
# 可选：输出平均结构
average crdset avg_backbone.rst
 
run
quit

聚类，三条轨迹合并（根据口袋周围5埃的关键残基进行聚类）

Pymol中确定残基名称
# 1. 定义口袋：选择 TRP 残基周围 5 埃内的所有完整残基
# （请将 123 替换为 TRP 的实际残基编号）
select pocket, byres resn TRP and resi 123 around 5
 
# 2. 显示名称
label pocket and name CA, "%s-%s" % (resn, resi)
# 3. 打印
>iterate (pocket and name CA), print("Residue: %s %s" % (resn, resi))

Cpptraj
parm pro.prmtop
resinfo :1-300
run
quit

Cpptraj.in（85.24）

parm com-dry.prmtop
 
trajin run_cmd_1/cmd-1.nc 1 last 10
trajin run_cmd_1/cmd-2.nc 1 last 10
trajin run_cmd_1/cmd-3.nc 1 last 10
trajin run_cmd_1/cmd-4.nc 1 last 10
trajin run_cmd_2/cmd-1.nc 1 last 10
trajin run_cmd_2/cmd-2.nc 1 last 10
trajin run_cmd_2/cmd-3.nc 1 last 10
trajin run_cmd_2/cmd-4.nc 1 last 10
 
trajin run_cmd_3/cmd-1.nc 1 last 10
trajin run_cmd_3/cmd-2.nc 1 last 10
trajin run_cmd_3/cmd-3.nc 1 last 10
trajin run_cmd_3/cmd-4.nc 1 last 10
 
autoimage #PBC处理
cluster C1 \
    hieragglo epsilon 0.6 linkage epsilonplot epsilonplot.dat \
    rms :13,15,16,19,70,71,72,90,91,92,115,144,146,194,195,196,199,219,220,223,236 \
    sieve 10 random \
    out cnumvtime.dat \
    sil Sil \
    summary summary.dat \
    info info.dat \
        cpopvtime cpopvtime.agr normframe \
    repout rep repfmt pdb \
    singlerepout singlerep.nc singlerepfmt netcdf \
    avgout Avg avgfmt restart
run
EOF

MM/GBSA[浩伦46] （85.24）

Input file for running PB and GB
&general
  startframe=0, endframe=100000, interval=200, keep_files=0,
/
&gb
  igb=5, saltcon=0.150,
/
&decomp
   idecomp=1, csv_format=0,
/
 
nohup mpirun -np 16 MMPBSA.py.MPI -O -i mmgbsa.in -o mmgbsa.dat -cp native.prmtop -rp pro.prmtop -lp lig.prmtop -y dry-align.xtc -eo mmgbsa.csv &
 
 
 
 
molecular dynamics "你的关键词" source:Journal source:of source:the source:American source:Chemical source:Society
在google scholar里用这个关键词搜
我一般看jacs，你也可以看看jctc，jmc，或者子刊的methods
 
 
 
 [hX1]GaMD在老版本amber中不支持
 [hX2]# 1. 检查 PDB 文件中是否记录缺失残基 
grep "^REMARK 465" input.pdb 
# 2. 检查 PDB 文件中是否记录缺失原子 grep "^REMARK 470" input.pdb 
# 3. 提取蛋白中实际存在的残基编号，推荐使用固定列宽方式
awk ' substr($0,1,4)=="ATOM"{ chain=substr($0,22,1) resi=substr($0,23,4) resn=substr($0,18,3) gsub(/ /,"",resi) key=chain" "resi if(key!=prev_key){ print chain, resi, resn prev_key=key } }' input.pdb > residue_list.txt
 # 4. 检查同一条链内是否存在残基编号跳跃 awk ' NR>1{ if($1==prev_chain && $2!=prev_res+1){ print "Possible gap:", prev_chain, prev_res, "->", $1, $2 } } { prev_chain=$1 prev_res=$2 }' residue_list.txt
 
 [LW3]看一下这个Se是否有特殊作用，特别是是否影响配体结合，如果有需要跟合作课题组讨论一下
 [浩伦4]第一，231 位是否在配体口袋附近。

 [浩伦38]Gaussian 16:  ES64L-G16RevB.01 20-Dec-2017