#Lisa Hall 5/6/10: calculates molecular g(r) #not including intramolecular correlations #thanks to Jennie Thomas for sharing her gofr script with vmd-l #so I had somewhere to start #loops over the different types of atoms in the molecules #for now I just have charge =0 and charge =-1 as the two types for {set t1 0} {$t1 < 2} {inc t1} { for {set t2 $t1} {$t2 < 2} {inc t2} { set outfile [open gofr$t1$t2.dat w] puts $outfile "r g$t1$t2" #initialize r and grsum by doing measure gofr from the first #molecule, "residue 0", to the others set sel1 [atomselect top "residue 0 and charge = -$t1"] set sel2 [atomselect top "residue != 0 and charge = -$t2"] set g [measure gofr $sel1 $sel2 delta .02 rmax 12 usepbc 1 first 1 last -1 step 1] set r [lindex $g 0] set grsum [lindex $g 1] #now loop over all other molecules, adding each g(r) to the total; #I have 100 molecules for {set i 1} {$i < 100} {incr i} { #find g_1(molecule i)-2(other molecules) set sel1 [atomselect top "residue $i and charge = -$t1"] set sel2 [atomselect top "residue != $i and charge = -$t2"] set g [measure gofr $sel1 $sel2 delta .02 rmax 12 usepbc 1 first 1 last -1 step 1] set gcurrent [lindex $g 1] set j 0 foreach gs $grsum gc $gcurrent { lset grsum $j [expr {$gs + $gc}] incr j } } set j 0 foreach gs $grsum { #renormalize: multiply by (N-1)/N and divide by N #since you have summed N g(r) #which were each normalized to 1/(N-1) lset grsum $j [expr {$gs*0.99/100.}] incr j } foreach jj $r kk $grsum { puts $outfile "$jj $kk" } close $outfile } }