Jump to content

Talk:Mod:Hunt Research Group/ESP charges for VMD

From ChemWiki

visualise charges using VMD

  • this allows us to use other colour scales by employing VMD to visualise the charges

set-up stuff

  • run extract_structure.py to form file.xyz
  • python path/extract_structure.py logfile (without .log extension)
  • save the script below into a file "charges_pdb.py"
  • this is a python script to combine standard orientation and charges into a pdb file with structure and charges in the correct columns for visualisation in VMD

run the script

  • run extract_charges.py to form file.chg
  • this will extract NBO and CHelpG charges
  • then run the charges_pdb.py script
  • run by typing "python charges_pdb.py arg1 file_name (without extension)
  • arg1 options are nbo or esp or both
  • nbo or both produces nbo charge file file_nbo.pdb
  • esp or both produces a chelpg charge file file_esp.pdb
  • note all files must have the same base name file_name

using VMD

  • read the pdb file into VMD then select as colouring method "occupancy"
  • best is to use a script to load the structure, the script load.tcl is given below
    • copy the pdb file you want to load to mol.pdb
    • in VMD type source load.pdb to activate the script in VMD
    • in VME then type "load" to run the script
    • this will generate a cpk molecule in the window
  • now manipulate the representation
    • click on tabs "graphics" and "representation"
    • in the representation window select "occupation" under colouring method
  • change the background colour
    • click on tabs "graphics" and "colours"
    • change the background colour if you want: Display, background
  • make sure the scale is where you want it
    • look in the pdb file, what is the absolute max/min, for example 0.29
    • click on tab "Trajectory"
    • change the "color scale data range" to the range you want
    • click set and watch the colors change

•you can make changes via the tcl console as well

    • mol scaleminmax 0 molID min max
    • use "molinfo list** to get information on the molecules
    • mol scaleminmax 0 0 -1.0 1.0
  • to improve the rendering use the following command
    • render Tachyon my_file "/Applications/tachyon/compile/macosx-x86-thr/tachyon" -aasamples 12 %s -mediumshade -res 1000 1000 -format TARGA -o %s.tga

charges_pdb.py

#
# charges_pdb.py       
#
# python script to combine standard orientation and charges into a pdb file  
# with structure and charges in the correct column for visualisation in VMD
# first run extract_structure.py to form file.xyz
# then run extract_charges.py to form file.chg 
# then run this script 
# nbo produces nbo charge file, esp produces chelpg charge file
# type python charges_pdb nbo/esp/both file_name (without extension)
# file_nbo.pdb or file_esp.pdb or both will be produced
# note all files must have the same base name
# read the pdb file into VMD then select as colouring method "occupancy"
#
import sys
import os
dir=os.getcwd()
#
if len(sys.argv) < 3:
  print('to run the script please type: \
     python charges_pdb.py nbo/esp/both input_file_name \(without extension\)')
  sys.exit()
else:
  base=str(sys.argv[2])
  file=base+'.log'
  geom=base+'.xyz'
  charge=base+'.chg'
  s='directory is '
  print('{0:}{1:}'.format(s,dir))
  s='xyz file is '
  print('{0:}{1:}'.format(s,geom))
  s='charges file is '
  print('{0:}{1:}'.format(s,charge))
#
  if sys.argv[1] == 'nbo' or sys.argv[1] == 'both' : 
    out1=base+'_nbo.pdb'
    s='output file is '
    print('{0:}{1:}'.format(s,out1))
  if sys.argv[1] == 'esp' or sys.argv[1] == 'both' : 
    out2=base+'_esp.pdb'
    s='output file is '
    print('{0:}{1:}'.format(s,out2))
  if sys.argv[1] != 'esp' and sys.argv[1] != 'nbo' and sys.argv[1] != 'both':
    print('to run the script please type: \
     python charges_pdb.py nbo/esp/both input_file_name \(without extension\)')
# endif
    
#close if
#
# setup 
total=0
col_num=[]
col_x=[]
col_y=[]
col_z=[]
col_atom=[]
col_esp=[]
col_nbo=[]

#==get structure==
f=open(geom,"r")
line =f.readline()
a1,b1,c1,d1=line.rstrip().split()
total=int(d1)
line =f.readline()
n=0
while n < total:
  line =f.readline()
  a1,b1,c1,d1=line.rstrip().split() 
  col_num.append(int(a1))
  col_x.append(float(b1))
  col_y.append(float(c1))
  col_z.append(float(d1))
  n=n+1
#close while
f.close()

#==get charges==   
f=open(charge,"r")
line =f.readline()
line =f.readline()
n=0
while n < total:
  line =f.readline()
  a1,b1,c1,d1=line.rstrip().split() 
  col_atom.append(b1)
  col_esp.append(float(c1))
  col_nbo.append(float(d1))
  n=n+1
#close while

#==write to .pdb file==
if sys.argv[1] == 'nbo' or sys.argv[1] == 'both' : 
  c=open(out1,"w")
  s='TITLE'
  c.write('{0:} \n'.format(s))
  s='REMARK'
  c.write('{0:} \n'.format(s))
  n=0
  s='ATOM'
  while n < total:
    c.write('{0:4}  {1:>5d}  {2:>3s}             {3:7.3f} {4:7.3f} {5:7.3f} {6:6.2f} \n'.\
    format(s,col_num[n],col_atom[n],col_x[n],col_y[n],col_z[n],col_nbo[n]))
    n=n+1
# close while
  s='END'
  c.write('{0:} \n'.format(s))
  c.close()
#endif

if sys.argv[1] == 'esp' or sys.argv[1] == 'both' : 
  c=open(out2,"w")
  s='TITLE'
  c.write('{0:} \n'.format(s))
  s='REMARK'
  c.write('{0:} \n'.format(s))
  n=0
  s='ATOM'
  while n < total:
    c.write('{0:4}  {1:>5d}  {2:>3s}             {3:7.3f} {4:7.3f} {5:7.3f} {6:6.2f} \n'.\
    format(s,col_num[n],col_atom[n],col_x[n],col_y[n],col_z[n],col_esp[n]))
    n=n+1
  #close while
  s='END'
  c.write('{0:} \n'.format(s))
  c.close()
#endif
#

load.tcl

proc load {args} {
#
#setup some display settings first
# set the background colour
# all atoms same size=orthographic depth=perspective
# turn off the depth shading
# turn off the axes
color Display Background blue2
display projection   Orthographic
display nearclip set 0.000000
display farclip  set 10.000000
display depthcue   off
axes location off
#
# read in the den cube file
mol new mol.pdb type pdb 
#
# set the molecule representation
# this is rep0
#
mol delrep 0 top
mol selection {all}
mol representation CPK 1.100000 0.300000 20.000000 16.000000
mol color Name
mol addrep top
#
# first reverse the colour scale so that red is positive and blue is negative
# mol modcolor rep_number mol_number coloring_method
# mol scaleminmax mol_number rep_number min max          
# note that the above mol vs rep are reversed!
# you will need to update the molecule_number as you load files!
color scale method BWR
mol modcolor 1 0 Volume 1 
mol scaleminmax 0 1 -0.050000 0.0500000
#
# you can add a colour scale bar from the extension menu manually
#
}
Retrieved from "https://chemwiki.ch.ic.ac.uk/index.php?title=Talk:Mod:Hunt_Research_Group/ESP_charges_for_VMD&oldid=734191"