c main program c c***purpose c c Computer simulation of supercoiled DNA by BD with analysis c of site juxtaposition process c c***description c c double precision version c c Program units: c ============= c length -- 10 nm, temperature -- K, concentration -- M, c energy -- kT, mass -- 1 for each bead. c c "common" blocks: c =============== c bf_coords axes of body-fixed coordiantes c ddx bond vectors c debyehu electrostatics c diffu_tensor hydrodynamic diffusion tensor c euler Euler angles c factor some constants c figinfo information of current configuration c ischeme DNA type c juxta_distribution information of juxtaposition process c jux_info information of juxtaposition process c matrix coords, bond angles c matrix2 mutual distances c myrandom BD variance c param simulation parameters c Pdl juxtaposition percentage c tor_param simulation parameters for cicular close DNA c xoldtw old configuration for checking c c Array g stores the bending rigidity for different subsegment c number up to 20. Typically, a choice of k >= 10 is good. c c c Input files c ========== c-------------------------------------------------------------------- c c a) File "Init_conf" c c The starting configuration. c For different 'initial' value (set in 'Input_data') c initial=1 -- no "Init_conf" is needed. c =2 -- "Init_conf" contains x(1:n,3) [coordinates] c =3/4 -- "Init_conf" contains x(),f(),u(),v() c [coordinates with local coordiante systems] c c-------------------------------------------------------------------- c c b) File "Init_fort" c c For coninued trajectory with juxtaposed information c Only needed when 'initial=4' in "Input_data". c c --Used to store number of juxtaposed for a certain separation length c in the trajectory to be continued; c --May need fort.7xx as supplemental input of juxtaposed information c --Check cont_juxta() at juxta.f for more details c c-------------------------------------------------------------------- c c c) File "Init_juxtadistance" c c This file is only needed when 'np600=1' in 'Input_data', for c site juxtaposition dissociation time output. c c njudge -- number of different criteria for dissociation c time c sjuxta_judge(1..njudge) -- value of those criteria c (unit: 10 nm) c c For more details, check init_juxta.f and out_dissociation_time.f c for the whole procedure recording and outputing dissociation time. c c-------------------------------------------------------------------- c c d) File "Init_juxtapair" c c This file is only needed when 'np600>=1' in 'Input_data'. c c For 'np600=1', output the pair-distance (Init_juxtapair) c This file should contain: c c1. PAIRS BEING FOLLOWED: c npair -- number of pairs to be explicitly followed c npair(1:npair,1:2) -- beads id for those pairs c c For 'np600=2' output the distribution about pair-distance c This file should contain: c c1. PAIRS BEING FOLLOWED: c npair -- number of pairs to be explicitly followed c npair(1:npair,1:2) -- beads id for those pairs c c2. DISSPCIATION CRITERIA: c njudge -- number of different criteria for dissociation c time c sjuxta_judge(1..njudge) -- value of those criteria c (unit: 10 nm) c c njudge c sjuxta_judge(1:njudge) c c3.BUCKETS FOR DISTRIBUTION OF DISTANCE BETWEEN TWO SITES: c d_min,d_max -- maximum and minimum values of distribution c -- to be outputed c (change nbucket in init_juxta.f for total number of buckets) c c For more details, check init_juxta.f, out_juxtapair.f(for c np600=1), and out_juxta_distribution(for np600=2). c c-------------------------------------------------------------------- c c e) File "Input_data" -- Most parameters needed for BD simulation c c1.SIZE c c n - number of beads c nc - subcharges (recommend 5 for 0.1M and subseg=10) c c2. TRAJECTORY: c c initial 1--initial from a circle c 2--read from 'Init_conf' without euler c 3--read from 'Init_conf' with euler data c 4--read from 'Init_conf' with euler data and juxtaposed c information c nbegin 1--from the first conformation c else--continue the trajectories with fort.117 c nfinal the last step number (nbegin..nfinal) c dtime time step (in ps) (recommmend 450 ps) c c3. SHAPE: c c icircular 1--circular c 0--linear C iclose 1--closed DNA (torsional force) c 0--nicked c iknot if knotting properties need to be checked c 1--check 0--no check c iele 0--no electrostatic interaction c 1--electrostatic interaction c c4. ALGORITHM: c ihydro 1--hydrodynamics c 0--no hydrodynamics c ibd 2--2 order BD c 1--1st order BD c note: hydrodynamics only have 2nd order version c nupdate the update frequency of hydrodynamics tensor c (recommend 10) c c5. RANDOM SEEDS: c rand1,rand2 c two integers between 1 and 2147483398 c c6. CONSTANT: c temp temparature (in K) (recommmend 298.0) c d dielectrostatic constant of solution c recommend (78.0 or 80.0) c visco viscosity (1.0 c statbl length of statistical length (10) c ss1 strentching energy (recommmend 100.0 ~100kT) c cg torsional constant (3.00) c phi0 natural torsion of DNA (10.5 bp/turn) c c 7. PARAMETERS: c c subseg segments in a statistical length (10~20) c bendc bending constant (correlated with subseg) c (subseg=10 bendc=4.806 subseg=20 bendc=9.820) c beadr hydrodynamic radius c (0.178 while subseg=20) c concen concentration of monoions (0-0.2M) c charge (e/10nm) c concen 0.01 0.02 0.05 0.1 0.2 0.5 1.0 c charge 24.3 29.6 41.5 60.8 99.4 262 788 c superdensity 0 -- -0.06 c c8. EXCLUDED VOLUME: c c interval exclude the electro interaction from i and i+intel-1 c segments (inteval >=1) c actually using (intel=inteval/r0) c ddna drmin parameters for artificial potential 1 c usually set ddna=0.22 (the range of wall) c drmin is artificial pushing (0.05~0.25) c c9. OUTPUT: c c np108,n108 --np108=1 output writhe number every c nwrithe steps (gyration) c np129,n129 --np129=1 output config information every c nprint steps c npbranch,nbranch--output the tip positions c c np600 --if pair-distance is interested c np600<0 nothing is interested c np600=0 output the dissociation time c np600=1 output the pair-distance (Init_juxtapair) c np600=2 output the distribution about pair-distance c neuler -- output (f,v,u) or not (recommend 1) c njuxta -- 1 output juxtaposition time c -- 0 no need to follow site juxtaposition c c-------------------------------------------------------------------- c c f) File "Input_juxta" c c This file is needed only when site juxtaposition time is calculated. c (njuxta=1 in "Input_data") c c1. SITE JUXTAPOSITION CRITERIA c c ndis - number of different criteria checked in a single trajectory c djux() - value of site juxtaposition criteria c should input exactly ndis values c (unit: 10 nm) c c2. SEPERATION LENGTH BETWEEN TWO SITES c c njux - number of different separation lengths checked in this c BD trajectory c iset() - value of separation lengths c (unit: beads) c c-------------------------------------------------------------------- c c***end implicit double precision(a-h, o-z) parameter(maxa=500,maxn=3*maxa,ncmax=5, maxb=10) parameter(pi=3.1415926535d0, ee=1.0d-5) integer num,nfinal,na,nb,nhb,igo,i,j double precision superdens common/debyehu/econst,debye,f_art,ddna,rmax,rmax2,rmin, +nc,nc3 common/ischeme/icircular,iele,ihydro,iclose,iknot,neuler common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma, +amass,am(maxn),coef,dtime2 common/tor_param/dr,cg,phi0,dlk,phii,subseg common/xoldtw/olddlk,xold(maxa,3),fold(maxa,3),vold(maxa,3) +,uold(maxa,3) common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3) common/matrix/x(maxa,3),bangle(maxa),rz(maxa) common/matrix2/ r(maxa,maxa) common/factor/dc1,dc2,bead2r,beadr2r,beadf1,beadf2 c****************************************************************** c INITIALIZATION c 1. read from Input_data open(unit=333,file='Input_data',status='old') read(333,*) n, nc, + initial, nbegin, nfinal, dtime, + icircular, iclose, iknot, iele, + ihydro, ibd, nupdate, + rand1, rand2, + temp, d, visco, statbl, ss1, cg, phi0, + subseg,bendc,beadr,concen,charge,superdensity, + interval,ddna,drmin, + np108,n108,np129,n129,npbranch,nbranch, + np600,neuler,njuxta close(333) c 2. constants setup r0=statbl/subseg ss1=ss1/(r0**2) econst=(charge*r0*48.0/nc)**2/(d*1.38*temp) dlk=superdensity*n*statbl*30.0/subseg/phi0 olddlk=dlk debye=dsqrt(3.246d+5*concen/temp) bpleng=0.034d0 bpmass=660.d0 fmass=statbl*bpmass/(bpleng*subseg) ftime=dsqrt(1.66*fmass/(1.38*temp)) dtime2=dtime dtime=dtime2*1.d-2/ftime visco=visco*ftime*1d+6/(fmass*1.66) cg = cg/(1.38d-3*temp*r0) c hydrodynamics factor dc2 = 1.d0/(8.d0*pi*visco) bead2r = 2.d0*beadr beadr2r = beadr*bead2r beadf1 = 9.d0/(32.d0*beadr) beadf2 = 3.d0/(32.d0*beadr) aa=beadr api=6*pi*visco*aa dc1=1.d0/api coef=dtime/api c artificial potential 1 f_art=drmin/coef c electrostatic trucation rmax=-log(ee*r0/econst/n/n)/debye rmax2=rmax+2*r0 nc3=3*nc c--------------------------------------------------------- c 3. other initializaion subroutines call init_juxta(n,np600) call init_rand(n,api) call ranset(nrand1,nrand2) call init_phi0(n) if (ihydro.eq.1) then call init_diffu(n) endif if (njuxta.eq.1) then call set_juxta(n) if (initial.eq.4) then initial=3 call cont_juxta(n) endif endif call setup_conf(n,initial) if (((icircular.eq.1).and.(iclose.eq.1)).and.(neuler.eq.1)) then call out_conf_e(n,0) else call out_conf(n,0) endif c set /xold/ do i=1,n do j=1,3 xold(i,j)=x(i,j) fold(i,j)=f(i,j) vold(i,j)=v(i,j) uold(i,j)=u(i,j) enddo enddo call update_rij(n) c****************************************************************** c SIMULATION PART ntime=nbegin do 1000 ntime0=nbegin, nfinal ntime=ntime+1 c--------------------------------------------------------- c1. BD propagation c c ihydro=0 no hydrodynamic c ibd=1 1st BD if ((ihydro.eq.0).and.(ibd.eq.1)) then call bd_nohydro1(n,ntime) c ibd=1 2nd BD elseif ((ihydro.eq.0).and.(ibd.eq.2)) then call bd_nohydro2(n,ntime) c ihydro=1 w. hydrodynamics always 2nd BD elseif (ihydro.eq.1) then if (mod(ntime,nupdate).eq.1) then call bd_hydro(n,ntime) else call bd_hy0(n,ntime) endif else print *, "Wrong selection of BD (with hydro). Stop." stop endif c--------------------------------------------------------- c 2. Output some properties if (np600.eq.0) then c call out_juxta_distance(n,ntime) call out_dissociation_time(ntime,n) endif if (np600.eq.1) then call out_juxtapair(ntime) endif if (np600.eq.2) then call out_juxta_distribution(ntime) endif if (njuxta.eq.1) then call juxta(n,ntime) endif if ((np108.eq.1).and.(mod(ntime, n108).eq.0)) then call cal_writhe(n,ntime,n108) endif if ((np129.eq.1).and.(mod(ntime,n129).eq.0)) then if (neuler.eq.1) then call out_conf_e(n,ntime) else call out_conf(n,ntime) endif c for debug, output the distribution of bond length smean=0.0;svar=0.0 do i=1,n-1 smean=smean+r(i,i+1) enddo smean=smean/(n-1) do i=1,n-1 svar=(r(i,i+1)-smean)**2+svar enddo svar=sqrt(svar/(n-1)) write(919,*) ntime,smean,svar call flush(919) call flush(129) endif if ((npbranch.eq.1).and.(mod(ntime,nbranch).eq.0)) then call out_tip(n,ntime) endif 1000 enddo c******************************************************************c c END Output the final configuration and also fort.117 to continue a c trajectory call out_final(n) call out_f117(n) stop end