solvate using cgwat

remember here waters are coarse grain, now the option -b should be 5.0

gromacs : coarse grain

need to prepare water explicitly
download dssp file for protein at :
ftp://ftp.cmbi.kun.nl/pub/molbio/data/dssp/

-bash: make: command not found(MAC)

Well, you need to download Xcode and install it on your mac.

compiling tcl hacked on NERSC Franklin(CRAY-XT3)

./configure --enhable-gcc --disable-load --prefix=$dir

Problem :

./../unix/tclUnixInit.c
./../unix/tclUnixInit.c: In function 'TclpSetVariables':
./../unix/tclUnixInit.c:536: error: storage size of 'name' isn't known
./../unix/tclUnixInit.c:547: warning: implicit declaration of function 'uname'
./../unix/tclUnixInit.c:536: warning: unused variable 'name'
make: *** [tclUnixInit.o] Error 1

Solution:
Add -DNO_UNAME to your CFLAGS in the Makefile

cd $dir/lib
mv libtcl8.3.so libtcl8.3.a

3d plot error bar in matlab

http://code.izzid.com/2007/Aug/Matlab_3D_Plot_with_Errorbars/

function [h]=plot3d_errorbars(x, y, z, e)
% this matlab function plots 3d data using the plot3 function
% it adds vertical errorbars to each point symmetric around z
% I experimented a little with creating the standard horizontal hash
% tops the error bars in a 2d plot, but it creates a mess when you rotate
% the plot
%
% x = xaxis, y = yaxis, z = zaxis, e = error value

% create the standard 3d scatterplot
hold off;
h=plot3(x, y, z, '.k');

% looks better with large points
set(h, 'MarkerSize', 25);
hold on

% now draw the vertical errorbar for each point
for i=1:length(x)
xV = [x(i); x(i)];
yV = [y(i); y(i)];
zMin = z(i) + e(i);
zMax = z(i) - e(i);

zV = [zMin, zMax];
% draw vertical error bar
h=plot3(xV, yV, zV, '-k');
set(h, 'LineWidth', 2);
end

error

/usr/lib64/libc.a(strtoul.o)(.text+0x0): In function `strtoul':
../sysdeps/generic/strtol.c:110: multiple definition of `strtoul'
/u0/b/bingo/tcl/tcl-crayxt3/lib/libtcl8.3.a(strtoul.o)(.text+0x0): first defined here
/usr/bin/ld: Warning: size of symbol `strtoul' changed from 521 in /u0/b/bingo/tcl/tcl-crayxt3/lib/libtcl8.3.a(strtoul.o) to 92 in
/u0/b/bingo/tcl/tcl-crayxt3/lib/libtcl8.3.a(tclUnixPipe.o)(.text+0x104): In function `TclpCreateTempFile':
: warning: the use of `tmpnam' is dangerous, better use `mkstemp'
/u0/b/bingo/tcl/tcl-crayxt3/lib/libtcl8.3.a(tclUnixFile.o)(.text+0xab7): In function `TclpGetCwd':
: warning: the `getwd' function is dangerous and should not be used.
child process exit status 1: /usr/bin/ld
Fatal Error by charmc in directory /u0/b/bingo/DEISA_BENCH/applications/NAMD/tmp/namd_Cray-XT4-HECToR_apoa1_i000006/n32p2t1_t001_i01/src/NAMD_2.6_Source/CRAY-XT3

Tcl ON Mac

Never try to screw with your tcl again. If already did, go to the following web, download the package and install it.(Don't try to compile tcl from scratch!!!!)

http://tcltkaqua.sourceforge.net/

remote access:
system preferece/sharing

interpolate missing residues in pdb file

1: download modeller(http://salilab.org/modeller/wiki/Missing%20residues)
2: you can get FASTA file from RCSB PDB(sequence details)
3: use sequence.py to generate .ali files
#################################
import sys
from modeller import *
from modeller.scripts import complete_pdb

# Get the sequence of the 1qg8 PDB file, and write to an alignment file
chain_id = sys.argv[1]

env = environ()
env.io.atom_files_directory = ['.', '../atom_files']
env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

aln = alignment(env)
mdl = model(env)

# read sequence from pdb list
code = 'XXXX'
mdl.read(file=code, model_segment=('FIRST:' + chain_id, 'LAST:' + chain_id))

# add the sequence to the alignment
aln.append_model(mdl, align_codes=code, atom_files=code)

# read complete pdb with missing residue
fasta='XXXX_' + chain_id + '.fasta.txt'
aln.append(file=fasta, alignment_format='FASTA')

# align them by sequence
aln.malign(gap_penalties_1d=(-500, -300))
aln.write(file='chain'+chain_id+'.ali')
#################################

4: use patch.py to add missing residues.

#################################
# insert the missing residues

import sys
from modeller import *
from modeller.automodel import * # Load the automodel class

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = ['.', '../atom_files']

chain_id = sys.argv[1]


class MyModel(automodel):
def select_atoms(self):
if chain_id == 'A' :
return selection(self.residue_range('817', '871'))
elif chain_id == 'C' :
return selection(self.residue_range('1', '7'),
self.residue_range('42', '61'),
self.residue_range('424', '431'))
elif chain_id == 'D' :
return selection(self.residue_range('1', '9'))
elif chain_id == 'E' :
return selection(self.residue_range('1', '8'),
self.residue_range('74', '76'))

a = MyModel(env, alnfile = 'chain' + chain_id + '.ali',
knowns = 'chain' + chain_id, sequence = 'chain'+chain_id+'_fill')
a.starting_model= 1
a.ending_model = 1

a.make()
#################################

periodic table

necessary components in protein translocation

mRNA: (generated in nucleus, transported to cytosol)
Another difference between eukaryotes and prokaryotes is mRNA transport. Because eukaryotic transcription and translation is compartmentally separated, eukaryotic mRNAs must be exported from the nucleus to the cytoplasm. Mature mRNAs are recognized by their processed modifications and then exported through the nuclear pore.

HOW THE PHOTOCOPIER WORKS

Before the advent of the photocopying machine, simple and cheap copying methods were both messy and crude. Older readers may remember the "Hectograph" from public school days. The master copy was written on a piece of paper with a special ink. The master was pressed face down on a tray of gelatine which dissolved and absorbed the ink. Moist blank sheets of paper were pressed onto the inked gelatine and a few inferior copies could be transferred this way.

The "Ditto" process is a variant on the Hectograph. By typing or writing, a special ink impression is made on the back of the master sheet. This master is mounted on a drum in a machine like a small printing press. On each rotation of the drum, methyl alcohol dissolves a bit of the ink and it is transferred to a blank sheet of paper. Twenty-five or so copies may be made.

The "Gestetner" process produces a master by typing, or writing with a stylus, on a special wax-impregnated paper. The wax is removed where the typewriter strikes it, making the paper porous to ink. Again, a small, rotary printing press is used. It is messy but can produce many copies.

In 1934, Chester Carlson, a physicist turned patent attorney, decided there must be a better way and embarked on a programme of research which led to the development of Xerography and the modern photocopier. The first commercial unit was the Zerox-914 which appeared in 1960.

There are a few substances such as selenium, arsenic, and tellurium which are "photoconductors." These substances are electrical insulators in the dark and are conductors in the light. A thin layer of this material is coated on a metal backing either in the form of a drum or a flat plate. At the start of the copy cycle, this drum is given a uniform static electric charge either by rubbing it or spraying on the charge.

Next, the page to be copied is illuminated with a strong light and an image of the page is formed by a lens on the charged photoconductor. White areas thus light up the photoconductor and it is rendered conducting so that its charge is leaked away to the metal backing. This leaves a latent image of the black areas of the page formed in the electric charge on the plate. This latent image must now be developed.

Originally, Carlson wafted a dust of black particles (called "toner") over the plate. Now small glass or plastic beads loaded with the black toner are either rolled or gently brushed over the plate. The toner sticks only where there is remaining charge; i.e., on the black areas of the original copy.

Next, a blank sheet of paper is rolled over the plate and the toner is transferred to the paper forming a perfect image of the original. This image, however, is very fragile and must be "fixed" or made permanent. This fixing is done in different ways, but one of the most common is by means of heat. The toner is actually a black coloured thermoplastic which, on being heated, melts and impregnates the paper. In other copiers, a chemical spray dissolves the toner so that it can impregnate the paper. Finally, the photoconductor plate is cleaned and ready to be recharged for the next cycle.

Machines in which the toner is wafted on, tend to deposit it at the edge of the charged regions so that large black areas are reproduced as black edges with white centres. The brushing method produces a more uniform toner deposition and thus reproduces large black areas more uniformly.

http://www.physics.uoguelph.ca/summer/scor/articles/scor54.htm

LDB problem

LDB: LOAD: AVG 4.34903 MAX 4.4356 MSGS: TOTAL 643 MAXC 13 MAXP 8 Refine

(http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l/2669.html)

Those LDB messages are part of every parallel run and are quite normal.
It's just reporting how the load balancer is doing. In this case, your
maximum load is only slightly more than your average load, and there are
171 total messages being sent between processors, which looks pretty good.

The older message from Brain Bennion was a hang during load balancing,
probably due to a compiler bug (a workaround is in 2.6b1).

Align molecule

(http://www.ks.uiuc.edu/Research/namd/tutorial/NCSA2002/hands-on/)
mol delete all

mol load psf pope.psf pdb pope_move.pdb
global popemol
set popemol [molinfo top]
mol load psf protein.psf pdb protein.pdb

display projection orthographic

mol modselect 0 $popemol "name P"
mol modstyle 0 $popemol VDW

puts {
     1. Switch to Move Molecule mode by moving the mouse into the graphics
        window and pressing the 8 key.  The cursor should now have the form
        of a cross-hair.

     2. While holding down the shift key, click on one of the protein atoms
        and drag the mouse.  You should see just the protein rotating
        around the atom you picked.  Try to align the axis of the protein
        with an imaginary line coming out of the screen, so that it's
        perpendicular to the membrane.  If you accidentally rotate the
        whole scene instead of just the protein, type "display resetview"
        in the text console to get back to the original view.

     3. Type "rotate x by -90" in the text console.  Now you want to center
        the protein between the phosphate atoms of the bilayer.  This time,
        click and hold on one of the protein atoms _without_ holding down
        shift.  Move the protein up and down until it's centered.

     4. Switch back to the orginal view with "display resetview".  If the
        protein looks good, you can go on to the next step; otherwise repeat
        steps 2 and 3 above.
}

Use pymol to build alpha helix

Build->Residue->Helix
Build->Residue->aa.

http://www.chem.ucsb.edu/~kalju/CSUPERB/public/pymol_movies.html

Click Mouse Mode to change from viewing to Editing

Easy to modify molecule using pymol(deleting)

**Nov 8**
In the Editing mode, Shift + mouse will move individual molecule.
Reset will move "disappearing molecules" to the screen.

making helix using Modeller

http://salilab.org/modeller/about_modeller.html

MODELLER is used for homology or comparative modeling of protein three-dimensional structures (1,2). The user provides an alignment of a sequence to be modeled with known related structures and MODELLER automatically calculates a model containing all non-hydrogen atoms. MODELLER implements comparative protein structure modeling by satisfaction of spatial restraints (3,4), and can perform many additional tasks, including de novo modeling of loops in protein structures, optimization of various models of protein structure with respect to a flexibly defined objective function, multiple alignment of protein sequences and/or structures, clustering, searching of sequence databases, comparison of protein structures, etc. MODELLER is available for download for most Unix/Linux systems, Windows, and Mac.

Using it, making helix becomes such a simple thing:
http://salilab.org/modeller/wiki/Make%20alpha%20helix
After installation, it's located at ///Library/modeller-9v5
terminal command: mod9v5

Color molecule using electrostatic potential

1: install apbs:
download the source code.
./configure --prefix=/use/local
make all
sudo make install

2: install pdb2pqr:(http://pdb2pqr.sourceforge.net/userguide.html#introduction)
(wasn't able to find numpy, so pdb2pka is disabled)
./configure --disable-pdb2pka --prefix=/usr/local
make
sudo make install

3: following the instructions on the tutorial
http://137.189.50.96/kbwong/teaching/pymol/pymol_tutorial.html

a. Use PDB2PQR to convert the PDB format to PQR format
> pdb2pqr.py --ff=amber --apbs-in 1w2i_nowat.pdb pymol.pqr
The PQR file will be output to pymol.pqr.

b. Use psize.py to determine the grid dimensions for APBS calculation
> psize.py pymol.pqr
You should be able to see the following results:

Center = 37.468 x 31.798 x 12.177 A
:
:
Coarse grid dims = 53.011 x 58.568 x 65.807 A
Fine grid dims = 51.183 x 54.452 x 58.710 A
Num. fine grid pts. = 97 x 97 x 97

Take a note on these parameters.

c. Edit the apbs.in
You need to enter the following parameters:

cgcent 37.468 31.798 12.177 # Grid Center
fgcent 37.468 31.798 12.177 # Grid Center

cglen 53.011 58.568 65.807 # coarse mesh lengths (A)
fglen 51.183 54.452 58.710 # fine mesh lengths (A)

dime 97 97 97 # Grid Points

d. Run APBS
> apbs apbs.in
After a while, it will create an electrostatics map called "pymol.dx".

just use vmd to import both .pdb and .dx file. Then choose drawing method to be surf, and coloring method to be volume. Then change the color scale data range in Trajectory bar to be like -8.00 to 8.00.

Beautiful plot appeared.!

blue surface : positively charged
red surface : negatively charged
white surface : neutral

Alpha Helix

http://www.web-books.com/MoBio/Free/Ch2C4.htm

An å helix has the following features:

• every 3.6 residues make one turn(100 degree, 1.5 A for each residue),
• the distance between two turns is 0.54 nm,
• the C=O (or N-H) of one turn is hydrogen bonded to N-H (or C=O) of the neighboring turn.

ß sheet:


Propensity of amino acid to form a alpha helix(from : Science [0036-8075] BLABER yr:1993 vol:260 iss:5114 pg:1637 -1640)

stop typing pwd using ssh

To use public keys with an ssh server, you'll first need to generate a public/private key pair:

$ ssh-keygen -t rsa
        

You can also use -t dsa for DSA keys, or -t rsa1 if you're using Protocol v1. (And shame on you if you are! Upgrade to v2 as soon as you can!)

After you enter the above command, you should see something like this:

Generating public/private rsa key pair.
Enter file in which to save the key (/home/rob/.ssh/id_rsa):

Just hit Enter there. It will then ask you for a pass phrase; just hit enter twice (but read the Security note below). Here's what the results should look like:

Enter passphrase (empty for no passphrase):
Enter same passphrase again:
Your identification has been saved in /home/rob/.ssh/id_rsa.
Your public key has been saved in /home/rob/.ssh/id_rsa.pub.
The key fingerprint is:
a6:5c:c3:eb:18:94:0b:06:a1:a6:29:58:fa:80:0a:bc rob@localhost

This created two files, ~/.ssh/id_rsa and ~/.ssh/id_rsa.pub. To use this keypair on a server, try this:

$ ssh server "mkdir .ssh; chmod 0700 .ssh"
$ scp .ssh/id_rsa.pub server:.ssh/authorized_keys2

easy way of compiling code?

http://www.deisa.eu/science/benchmarking

wget \
http://www.deisa.eu/science/benchmarking/files/benchmark_NAMD.tar.gz
gunzip benchmark_NAMD.tar.gz
tar xf benchmark_NAMD.tar
#--- Download NAMD source from http://www.ks.uiuc.edu/Research/namd/
mv NAMD_2.6_Source.tar.gz DEISA_BENCH/applications/NAMD/src/
cd DEISA_BENCH/applications/NAMD/
perl ../../bench/jube -debug bench-Cray-XT4-HECToR.xml

Protein terminus

Due to the chemical structure of the individual amino acids, the protein chain has directionality. The end of the protein with a free carboxyl group is known as the C-terminus or carboxy terminus, whereas the end with a free amino group is known as the N-terminus or amino terminus.

Lipid Tail Selection

(From dechandl@ks.uiuc.edu)
"
I just wanted to point out that selections like "lipids" in VMD are really
just pre-defined atomselect macros. You can create your own macros for more complex selections like this. For example, here is my macro for dealing with this for POPE, POPG and POPC - you'll have to change the selections if you have different lipid preferences:

atomselect macro lipidtails {(resname POPE and not name N HN1 HN2 HN3 C12 H12A H12B C11 H11A H11B P O11 O12 O13 O14 C1 HA HB C2 HS O21 O22 C21 C3 HX HY O31 O32
C31) or (resname POPC and not name N C12 C13 C14 H21 H22 H23 H31 H32 H33 H41 H42 H43 C11 H11 H12 C15 H51 H52 P1 O1 O2 O3 O4 C1 HA HB O21 O22 C21 C3 HX HY O31 O32 C31) or
(resname POPG and not name C12 O12 H31 H21 H22 C11 H11 O11 H12 C15 H51 H52 P1 O1 O2 O3 O4 C1 HA HB C2 HS O21 O22 C21 C3 HX HY O31 O32 C31)}

I keep these in my .vmdrc file so that they're defined at startup. Then selections like "lipidtails" will work from then on. This is useful if you know you'll want to re-use a selection often.
"

Take a look at the paper : bilayer.pdf(http://www.ks.uiuc.edu/Training/CaseStudies/pdfs/) to figure out how to choose the lipid head.

Make code working on NERSC

compile LAPACK on NERSC:
tar -zxvf lapack-3.0.tgz
cd LAPACK
cd INSTALL
cp make.inc.LINUX ../make.inc
cd ..
cd BLAS/SRC
make
cd ../..
make

VMD movie plugin for MAC

convert command comes from : ImageMagick
http://www.imagemagick.org/script/binary-releases.php
sudo port install imagemagick(from darwin port)

ppmtompeg comes from : Netpbm
http://netpbm.sourceforge.net/
sudo port install netpbm(from darwin port)

Wham calculation

Histogram must overlap somehow!!!! Otherwise, you would get really weird free energy value.

Hydrophobic effect

1:hydrophobic force:

A very frequently cited example of entropic force is hydrophobic force. Though hydrophobic force has an enthalpic component, about half of it originates from the entropy of the hydrogen bonded three dimensional network lattice of water molecules at room temperature. Since each water molecule is capable of donating two hydrogen bonds through the two protons and accept two more hydrogen bonds through the two sp3 hybridized lone pairs, a water molecule, unlike Hydrogen fluoride (accept 3 but donate only 1) or ammonia (donate 3 but accept only 1) which mainly form linear chains, can form an extended three dimensional network lattices. Introduction of a non-hydrogen bonding surface disrupts this networks and the rest of the network immediately rearranges itself around that surface so as to minimize the number of disrupted hydrogen bonds. If the introduced surface had ionic or polar nature there would have been hydrogen molecules standing more or less orthogonal to this surface. But a non-hydrogen bonding surface forces the surrounding hydrogen bonds to be tangential and become locked in a clathrate like basket shape. Water molecules involved in this clathrate like basket envelope formation around the non-hydrogen bonding surface are entropically restrained. Thus any event that would minimize this type of surface would be entropically favoured. For example when two such particles comes very close they merge and the water trapped in between would be released from the trapped state and join the free bulk water molecules and lead to increase in entropy. That is the basis of the so called "attraction" between hydrophobic objects in solution. Though when two hydrophobic objects attract when they are very close, there may be a mild repulsion when they are about to come close. This mild repulsion may happen because the water molecules when trapped in between two hydrophobic surfaces coming from different direction are frustrated about whether to orient tangential to one surface or the other. However when the surfaces come further close this repulsion is more than compensated by the great increase in entropy when the water molecules are completely displaced of the clathrate-like sheath.

/* =============================================== */
Chandler, D., "Insight Review: Interfaces and the driving force of hydrophobic assembly” Nature 437, 640-647 (2005).

Non-integer calculation in Shell

use bc for this:
example:

#!/bin/bash
# Linux Shell Scripting Tutorial 1.05r3, Summer-2002
# Written by Vivek G. Gite 
# Latest version can be found at http://www.nixcraft.com/
# Q10
a=5.66
b=8.67
c=`echo $a + $b | bc`
echo "$a + $b = $c"
# ./ch.sh: vivek-tech.com to nixcraft.com referance converted using this tool
# See the tool at http://www.nixcraft.com/uniqlinuxfeatures/tools/

Also, to run a loop of jobs, use the command "wait" to wait until the previous job finishes.( you need to run the job in background using "&")

problem with the solvation

I have the opposite z direction with Todd, so the box size in z direction should be inversed.

The deletion of lipid didn't work so well. So we should first minimize the system, to diminish the overlap of lipid and protein. Then from this structure, we generate the restrained structure.

wall -- send a message to everybody's terminal.

(From : http://linux.die.net/man)

Synopsis

wall [-n] [ message ]

Description

Wall sends a message to everybody logged in with their mesg(1) permission set to yes. The message can be given as an argument to wall, or it can be sent to wall's standard input. When using the standard input from a terminal, the message should be terminated with the EOF key (usually Control-D).

The length of the message is limited to 20 lines. For every invocation of wall a notification will be written to syslog, with facility LOG_USER and level LOG_INFO.

Options

-n

Suppresses the normal banner printed by wall, changing it to "Remote broadcast message". This option is only available for root if wall is installed set-group-id, and is used by rpc.walld(8).

Environment

Wall ignores the TZ variable - the time printed in the banner is based on the system's local time.

/* ========================================== */

write - send a message to another user

Synopsis

write user [ttyname]

Description

Write allows you to communicate with other users, by copying lines from your terminal to theirs.

When you run the write command, the user you are writing to gets a message of the form:

Message from yourname@yourhost on yourtty at hh:mm ...

Any further lines you enter will be copied to the specified user's terminal. If the other user wants to reply, they must run write as well.

When you are done, type an end-of-file or interrupt character. The other user will see the message EOF indicating that the conversation is over.

You can prevent people (other than the super-user) from writing to you with the mesg(1) command. Some commands, for example nroff(1) and pr(1), may disallow writing automatically, so that your output isn't overwritten.

If the user you want to write to is logged in on more than one terminal, you can specify which terminal to write to by specifying the terminal name as the second operand to the write command. Alternatively, you can let write select one of the terminals - it will pick the one with the shortest idle time. This is so that if the user is logged in at work and also dialed up from home, the message will go to the right place.

The traditional protocol for writing to someone is that the string '-o', either at the end of a line or on a line by itself, means that it's the other person's turn to talk. The string 'oo' means that the person believes the conversation to be over.

swig wrap c++ and tcl(Makefile)

TOP = .
SWIG = /home/bingo/bin/swig/bin/swig
CXXSRCS = example.cxx
TARGET = example
INTERFACE = example.i
LIBS = -lm -llapack (make sure to change the lib here, important)

all::
$(MAKE) -f $(TOP)/Makefile.in CXXSRCS='$(CXXSRCS)' SWIG='$(SWIG)' \
LIBS='$(LIBS)' TARGET='$(TARGET)' INTERFACE='$(INTERFACE)' tcl_cpp

static::
$(MAKE) -f $(TOP)/Makefile.in CXXSRCS='$(CXXSRCS)' SWIG='$(SWIG)' \
TARGET='mytclsh' INTERFACE='$(INTERFACE)' tclsh_cpp_static

clean::
$(MAKE) -f $(TOP)/Makefile.in tcl_clean

check: all

(if the compiling ends normally, but tcl doesn't run, it might means you are missing some flag during the make)

Call Lapack in C code:

use the grid class
http://oldmill.uchicago.edu/~wilder/Code/grid/

What is interesting about this C++ grid class?

• Simplicity. The class is (really) easy to build and use. This code compiles and runs fine without using any external libraries, but it can be compiled with optimized BLAS/LAPACK libraries to optimally run matrix computations.
• Speed. Many libraries are very clever in the way they optimize matrix calculations to avoid unnecessary copying and to carefully arrange the order of computations. My experience is that speed considerations are either irrelevant, in which case the optimizations are also not important, or critical, in which case other libraries often still do not get it quite right. This class provides simple interfaces to BLAS/LAPACK routines that can be called directly for these critical computations.
• Familiarity. The syntax is, in many ways, similar to what is used in 'R', and it is easy for people who are not experts in C++ to use.

OR, you can use the library(GSL):
http://www.gnu.org/software/gsl/

mutate protein in VMD

(From VMD mail list : http://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/2018.html)

This is a frequently requested feature. Here's the current solution.

If you're building a psf file using psfgen there is a mutate command that
can by applied during segment generation, for example:

segment BPTI {
pdb output/6PTI_protein.pdb
mutate 23 ALA ; # mutate residue 23 to ALA
}

When you try to load coordinates using coordpdb you'll get correct mapping
for the backbone atoms, but the sidechain may have some unfortunate name
conflicts. If this happens you can rename the atom names for that one
residue in the source pdb file (pdbalias won't help since it works on all
residues of a type). While you're at it, you can rename the residue so
you won't need the mutate command either. Here's a full example:

mol pdbload 6PTI
set res29 [atomselect top "resid 29"]
$res29 set resname ARG
set good [atomselect top "protein and not (resid 29 and not backbone)"]
$good writepdb /tmp/bpti.pdb
package require psfgen
topology /Projects/namd2/toppar/top_all22_prot.inp
segment BPTI {
pdb /tmp/bpti.pdb
}
coordpdb /tmp/bpti.pdb BPTI
guesscoord
writepsf bpti.psf
writepdb bpti.pdb
mol load psf bpti.psf pdb bpti.pdb

In this case we threw out all of the sidechain coordinates for the mutated
residue to be absolutely safe. This results in a different rotamer for
the side chain. If we use all of the coordinates (i.e., just "protein"
instead of "protein and not (resid 29 and not backbone)" we get a closer
match to the original structure since there is some uniformity in atom
names between different residue types. In a more complex case you could
use VMD's atom selections to rename individual atoms in the mutated
residue to match the corresponding atoms in the topology file.

matlab

%%% beautiful plot:
Linespec

plot(t,sin(2*t),'-mo',...
'LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor',[.49 1 .63],...
'MarkerSize',12)

%%%error bar
Matlab's plot command is not as nice as Gnuplot plot, and, for example, you can't use it to draw a plot with error bars. But drawing plots like that is so important that Matlab has a separate command for that. The command is errorbar , and you should use it in combination with plot, as shown below:

>> plot (x, y_prime)
>> hold on
>> errorbar(x, y, sigma_y, '.')
>> hold off
>>

The way it works is as follows. First we plot y' = a + bx versus xusing the plot command. Then we tell Matlab to hold on to the plot and add other elements to it rather than overwrite it with a new plot. Then we plot errorbars. The last argument in the errorbar command tells Matlab that the points given by (x_i, y_i) pairs should not be connected. Matlab draws them connected by default. The last command, hold off tells Matlab that the plot is finished, and no new elements will be added to it.

change vmd background color

Open Graphics->Colors, then select, Display->Background->White.

Line of Best Fit For Points in Three Dimensional Space

Theory:
(from : http://mathforum.org/library/drmath/view/69103.html)

The line you are looking for is called the 3D Orthogonal Distance Regression (ODR) line.  This line can be computed without writing a
lot of computer code if you already have code for the Eigen decomposition or the simgular value decomposition (SVD).

Let's parameterize the line like this,

x = x0 + at
y = y0 + bt
z = z0 + ct

where (x0,y0,z0) is a point on the line and (a,b,c) is a unit vector.

We want to minimize the sum of squared distances from the points to the line.  The vector equation for the sum leads to

f(x0,y0,z0,a,b,c) = SUM{[c(yi - y0) - b(z - z0)]^2 +
       [a(zi - z0) - c(x - x0)]^2 +
       [b(xi - x0) - a(y - y0)]^2}

If we take the first derivatives with respect to x0, y0, and z0 and set the results equal to zero we get equations that can be manipulated
to yield

(x0-xbar)/a = (y0-ybar)/b = (z0-zbar)/c

where (xbar,ybar,zbar) is the centroid of the data.  Notice that (xbar,ybar,zbar) is a valid choice for (x0,y0,z0).  This proves that
the centroid is on the ODR line, and so we choose it for (x0,y0,z0).

Now we just need to find vector (a,b,c).

We need a little notation to make the next steps easier.  Let C be the centroid, let L be the ODR line, and let P be the plane through C such
L is perpendicular to P.  We are doing this in order to make use of the ODR plane analysis, which is closely related.

By the Pythagorean theorem, for a set of points Xi,

SUM[distance(Xi,L)^2] = SUM[distance(Xi,C)^2] -
       SUM[distance(Xi,P)^2]

We want to choose (a,b,c) so that SUM[distance(Xi,L)^2] is minimized. This is the same as maximizing SUM[distance(Xi,P)^2].  (Notice that
SUM[distance(Xi,C)^2] is a constant).

The ODR plane analysis can be found here:

Orthogonal Distance Regression Planes
http://mathforum.org/library/drmath/view/63765.html

The ODR plane analysis uses matrices A and M.  For the ODR plane we want the eigen vector of A that corresponds to its smallest
eigenvalue.  This is the same as the singular vector of M that corresponds to its smallest singular value.  For the ODR line we
instead want the eigen vector of A that corresponds to its largest eigenvalue, or the singular vector of M that corresponds to its
largest singular value.  This maximizes the desired sum.  The decomposition of A or M gives us the vectors for both the ODR plane
and the ODR line at the same time.

To summarize, the 3D ODR line contains the centroid, and has as its
direction vector the eigen vector (or singular vector) referred to above

Shortest Distance between two Line Segments

http://homepage.univie.ac.at/Franz.Vesely/notes/hard_sticks/hst/hst.html

superimpose 2 plots together

Use Excel to import the picture.
Use gnuplot to get the data plot, save as a pdf file, and again, use Excel to import the pdf file. pdf is transparent in Excel!! So, now you can use mouse to adjust your plot and to overlap them together.

save the chart as a picture to print

gnuplot

LogScale:
set log x 2

change line width
plot "*.dat" w lp lw(linewidth) 10 ps(pointsize)

external links:
http://sparky.rice.edu/~hartigan/gnuplot.html

demo scripts:
http://gnuplot.sourceforge.net/demo/

http://t16web.lanl.gov/Kawano/gnuplot/index-e.html

/* ================================= */
%%% modify value while ploting

plot "accum_1histo_lg_12.0_pp_20.0.dat_bak" u ($1*0.0001):2 w l

write pdb file from dcd trajectory

mol load psf yourprotein.psf dcd yourprotein.dcd

set nf [molinfo top get numframes]

for { set i 0 } {$i < $nf } { incr i } {
set sel [atomselect top protein frame $i]
$sel writepdb $i.pdb
}

/* ============================== */
want to convert NAMD binary out to pdb file?
just wrong 0 step namd simulation. :-)

Units

pN : pico Newton.
NAMD force unit : 1 kcal/mol*A
1cal = 4.184J

NAMD

NAMD can use the beta coloumn of the PDB file to do the Constraints and Restraints calculation.
But the SMD use the beta column to specify which atom is fixed, and occupancy column for the SMD atoms.

SMD output
SMD prints out the timestep, current position of the Center of Mass of the restrained atoms and the current force applied to the center of mass (pN).