VMD-L Mailing List

From: Frank Zack (frankzack123_at_REMOVE_yahoo.de)
Date: Thu Nov 26 2015 - 06:30:57 CST

Hi all,
i've written a simple,dirty tcl-script to coarse-grain 40 Chains of EPR Poylmer (All-Atom,1000-repeating units each chain) to a Martini-like representation.It perfectly works for the first 4 chains...however during coarsegraining the 5th chain vmd gets killed and i have no idea why...
Do you have any suggestions?

proc cg {} {

        #package require pbctools
        #package require topotools 1.5 
        #topo readlammpsdata stage21.data

        set sel [atomselect 0 "residue 40 to 80"]
        set r1 [$sel get residue]
        set firstRes [lindex $r1 0]
        set lastRes [lindex $r1 end]
        set oldcell [pbc get]
        puts "firstRes:$firstRes"
        puts "lastRes:$lastRes"

        #always try to free mem to prevent vmd from getting killed
        $sel delete

        set counter     0

        set cname "1"
        set ppname  "13"

        #create new mol+atoms before           
        set sel_c  [atomselect 0 "name '$cname'" ]
        set sel_pp [atomselect 0 "name '$ppname'" ]
        set num_c  [$sel_c  num]
        set num_pp [$sel_pp num]
        set nrbeads [expr (($num_c/2)+$num_pp)]

        $sel_c  delete
        $sel_pp delete

        set mol [mol new atoms $nrbeads]
        animate dup $mol        set logfile [open "cgepr.log" "w"]
        for {set curchain $firstRes} {$curchain <= $lastRes} {incr curchain} {
                puts "coarsing chain $curchain"
                set sel [atomselect 0 "residue $curchain"]
                set maxid [lindex [$sel list] end]
                set maxid [expr ($maxid-1)]
                set bbatoms {$cname $ppname}
                set bb [atomselect 0 "residue $curchain and name '$cname' '$ppname'"]
                set bblist [$bb list]

                $bb delete
                set curid [lindex [$sel list] 0]
                $sel delete

                puts "maxid:$maxid"
                puts "curid:$curid"

                while {$curid < $maxid} {
                        #ADD BEAD
                        set cursel [atomselect 0 "index $curid"]
                        set nextid [expr ($curid+3)]
                        set nextsel [atomselect 0 "index $nextid"]
                        set curname [$cursel get name]
                        set nextname [$nextsel get name]
                        set cg  0

                        #ignore-terminators
                        set test [lsearch $bbatoms $curname]
                        if {$test==-1} {incr curid;continue}
                        if {$curname==$ppname} {incr curid;continue}

                        if {$curname==$cname} {
                                if {$nextname==$cname} {
                                        if { $nextid>$maxid } {break}
                                        set tmp [atomselect top "index $curid to $nextid"]
                                        set cg [measure center $tmp weight {12 1 1 12}]

                                        incr counter
                                        set sel2 [atomselect $mol "index [expr ($counter+1)]"]
                                        $sel2 set name "PE"
                                        $sel2 set type "1"
                                        $sel2 set {x y z} $cg
                                        $sel2 set resid $curchain

                                        incr curid
                                        incr curid
                                        incr curid
                                        $sel2 delete
                                        $tmp  delete
                                        continue

                                }
                                if {$nextname==$ppname} {
                                        if { $nextid>$maxid }
                                        #overnext ist das methyl am pp
                                        set overnext [expr ($nextid+1)]
                                        set tmp [atomselect 0 "index $curid to $overnext"]
                                        set cg [measure center $tmp weight {12 1 1 12 1 12}]
                                        #lappend namelist  "PP"
                                        #lappend typelist   2
                                        #set xyz [list $cg]
                                        #lappend xyzlist     $cg
                                        #lappend idlist      $counter
                                        #lappend reslist     $curchain
                                        #lappend resiudelist $curchain
                                        incr counter

                                        set sel2 [atomselect $mol "index $counter"]
                                        $sel2 set name "PP"
                                        $sel2 set type "2"
                                        $sel2 set {x y z} $cg
                                        $sel2 set resid $curchain

                                        incr curid
                                        incr curid
                                        incr curid
                                        $sel2 delete
                                        $tmp  delete
                                        continue
                                }
                        }
                }
        }
        #puts "adding $counter atoms"
        #set mol [mol new atoms $counter]
        #animate dup $mol

        #set sel [atomselect $mol all]
        #$sel set name $namelist
        #$sel set type $typelist
        #$sel set {x y z} $xyzlist
        #$sel set resid $reslist

        #ADDBONDS

        puts "adding bonds"
        puts "first:$firstRes, last:$lastRes"
        for {set i $first } {$i <=$last } {incr i} {
                #puts "i:$i"
                set sel [atomselect $mol "all resid $i" ]
                set l [$sel get index]
                set firstIndex [lindex $l 1]
                set lastIndex  [lindex $l end]
                puts $logfile "outerloop for chain $i"
                puts $logfile "firstIndex: $firstIndex"
                puts $logfile "lastIndex:  $lastIndex"
                for {set k $firstIndex} {$k <=$lastIndex} {incr k} {
                        set i1 $k
                        set i2 [expr ($k-1)]
                        ::TopoTools::addbond $mol $i1 $i2 1 1
                        puts $logfile "added bond for chain $i, index1:$i1, index2:$i2"
                }
        }
        set sel [atomselect $mol all]
        puts $logfile "guessing angles"
        ::TopoTools::guessangles $sel
        #pbc set cell [lindex $oldcell 0]
        ::PBCTools::pbcset [lindex $oldcell 0]
        topo writelammpsdata -sel $sel cg.data