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).