From: Carlo Guardiani (carlo.guardiani_at_dsf.unica.it)
Date: Tue Apr 24 2018 - 12:52:20 CDT
Dear NAMD experts,
I am trying to set up a simulation of a thin
layer of amorphous silicon dioxide with a hole
in the middle. Basically I am following the
tutorial "Modeling Nanopores for Sequencing
DNA" by Comer, Wells and Aksimentiev.
The dimensions of my system that I store in
file grid_basis.txt are:
49.78 0.0 0.0
0.0 49.78 0.0
0.0 0.0 6.948
0.0 0.0 0.0
Since I want to generate a pore parallel to the
z axis with diameter of 10 Ang I call the
gridSourcePore function as:
./grid/gridSourcePore grid_basis.txt 2 6.948 10 0 pore_points.txt
The output is the following:
Scanning points.
Total points: 2500
Source points: 64
Remaining points: 2436
Source fraction: 0.0256
Total volume: 17217.48028
Source volume: 440.7674952
Remaining volume: 16776.71279
The system is then built with the buildSystem.tcl
script that I customized like this:
-------------------------------------------------------------------------
set targetDensity 1.8
set remainingVolume 16776.71279
set nu [expr {$targetDensity*$remainingVolume/240.337}]
set n [expr {int(ceil(sqrt($nu)))}]
inorganicBuilder::initMaterials
set box [inorganicBuilder::newMaterialBox SiO2 {0 0 0} [list $n $n 1]]
inorganicBuilder::buildBox $box sio
--------------------------------------------------------------------------
I now call the distributeAtoms.tcl script to open
the pore:
--------------------------------------------------
mol delete all
# Parameters:
set poreLength 6.948
set poreDiameter 10.0
set poreAngle 0.0
set basisFile grid_basis.txt
set psf sio.psf
set pdb sio.pdb
set outPdb sio_ready.pdb
set pi [expr {4.0*atan(1.0)}]
set halfLen [expr {0.5*$poreLength}]
set s0 [expr {0.5*$poreDiameter}]
set slope [expr {tan($poreAngle*$pi/180.0)}]
proc inPoreWalls {pos} {
global halfLen s0 slope
foreach {x y z} $pos { break }
if {abs($z) > $halfLen} { return 0 }
set s [expr {sqrt($x*$x + $y*$y)}]
if {$s < $s0 + $slope*abs($z)} { return 0 }
return 1
}
# Load the basis vectors.
set cellVecList {}
set in [open $basisFile r]
foreach line [split [read $in] "\n"] {
lappend cellVecList [concat $line]
puts "cellBasisVector $line"
}
close $in
set cellA [lindex $cellVecList 0]
set cellB [lindex $cellVecList 1]
set cellC [lindex $cellVecList 2]
# Load the molecule.
mol load psf $psf pdb $pdb
set all [atomselect top all]
set num [$all num]
# Reposition the atoms.
puts "Repositioning atoms."
set posList {}
for {set i 0} {$i < $num} {incr i} {
set bad 1
# Find positions within the membrane.
while {$bad} {
set a [vecscale [expr {rand()-0.5}] $cellA]
set b [vecscale [expr {rand()-0.5}] $cellB]
set c [vecscale [expr {rand()-0.5}] $cellC]
set r [vecadd $a $b $c]
if {[inPoreWalls $r]} {
set bad 0
}
}
lappend posList $r
}
# Set the positions and write the pdb.
$all set {x y z} $posList
$all writepdb $outPdb
--------------------------------------------------------
Finally, I set charges and types for annealing:
----------------------------------------
mol load psf sio.psf pdb sio_ready.pdb
set all [atomselect top "all"]
set sil [atomselect top "type SI"]
$sil set charge 2.4
set oxy [atomselect top "type O"]
$oxy set charge -1.2
$oxy set type OSI
$all writepsf sio_ready.psf
#Mark the atoms for gridforce
$all set beta 1.0
$all set occupancy 1.0
$all writepdb sio_all.pdb
----------------------------------------
So far I didn't have any problem. However, when
I ran the minimization as
namd2 sio_anneal_min.namd > sio_anneal_min.log
the pore radius widened from 10 Ang to about
16 Ang. When I then ran the annealing
namd2 sio_anneal.namd > sio_anneal.log
my system was completely deformed from a slab
with a pore to something that looks like a oval
doughnut.
When I tried to re-run the simulations disabling the
gridforce, at the end of the minimization the system had
retained its diameter of 10 Ang but during the annealing
the pore closed and at the end I got an irregular disk
with no pore at all. Could you please tell me what to do ?
Is there a way to preserve the pore without forcing the
shape of the system to become hexagonal ? I would like to
retain the initial rectangular shape of my slab.
Many thanks for your help and kind regards,
Carlo Guardiani
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:21:02 CST