Skip to content

Add Reaction Dynamics to BD Code #19

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sarahsandholtz wants to merge 12 commits into
base: master
Choose a base branch
from
Open

Add Reaction Dynamics to BD Code #19

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
b932823
First version: tracking methylation profile during BD simulations
sarahsandholtz Aug 18, 2016
c2976c2
Merge remote-tracking branch 'upstream/master'
sarahsandholtz Aug 18, 2016
def90e8
fixed bugs: seg faults, i/o issues
Aug 20, 2016
fae0e26
fixed bugs: declaring variables and passing them between subroutines
sarahsandholtz Aug 30, 2016
8d828cc
added new variable and updated output to screen
sarahsandholtz Aug 30, 2016
01bc2ae
fixed bug in how methylation profile was updated
Aug 31, 2016
8b49ea8
fixed error in selecting site where methyl mark will spread
Sep 1, 2016
3b9d011
fixed bug with demethylation scheme
Sep 11, 2016
f2adedf
fixed bugs with modified timestep and demethylation scheme
Sep 21, 2016
42c4f10
fixed bug with finding which pairs of beads could react; also fixed p…
Oct 8, 2016
0fe9b6b
moved methyl profile code from Brownian Dynamics subroutine into main…
Oct 10, 2016
f84fb39
Moved methylation algorithm back into BD subroutine; added toggle to …
Oct 26, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
315 changes: 315 additions & 0 deletions BasicWLC/BDcode/BDsim.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,315 @@
!---------------------------------------------------------------*

SUBROUTINE BDsim(R,U,NT,N,NP,TIME,TTOT,DT,BROWN, &
INTON,IDUM,PARA,SIMTYPE,HAS_COLLIDED,FPT_DIST, &
COL_TYPE, METH_STATUS, KM, KD, NUM_SPREAD, IN_RXN_RAD, PAIRS, NUC_SITE, NUM_METHYLATED, NUM_DECAY)

!
! External subroutine to perform a Brownian dynamics simulation.
!
! Andrew Spakowitz
! Written 11-11-13

use mt19937, only : rnorm

DOUBLE PRECISION R(NT,3) ! Bead positions
DOUBLE PRECISION U(NT,3) ! Tangent vectors
DOUBLE PRECISION TIME ! Time of BD simulation
DOUBLE PRECISION TTOT ! Final time of BD simulation
INTEGER N,NP,NT ! Number of beads

! Variables in the simulation

DOUBLE PRECISION B(NT,1) ! Bond length
DOUBLE PRECISION RS(NT,3) ! R during the step
DOUBLE PRECISION US(NT,3) ! R during the step
DOUBLE PRECISION L0 ! Bond distances
DOUBLE PRECISION DT,DT0 ! Time step size
INTEGER RK ! Runge-Kutta index
DOUBLE PRECISION DRDT(NT,3,4) ! Position rate of change
DOUBLE PRECISION DUDT(NT,3,4) ! Position rate of change
INTEGER I,J,IB ! Index Holders
DOUBLE PRECISION DOTU
DOUBLE PRECISION R0(3)

! Variables for use in the force calculations

DOUBLE PRECISION FELAS(NT,3) ! Elastic force
DOUBLE PRECISION FPONP(NT,3) ! self-int force
DOUBLE PRECISION TELAS(NT,3) ! Elastic force
DOUBLE PRECISION TPONP(NT,3) ! self-int force
DOUBLE PRECISION FORCE ! External force
INTEGER FON ! Is force on?

! Variables in the simulation

DOUBLE PRECISION EB,EPAR,EPERP
DOUBLE PRECISION GAM,ETA
DOUBLE PRECISION XIR,XIU
DOUBLE PRECISION LBOX ! Box edge length
DOUBLE PRECISION LHC ! Length of HC int
DOUBLE PRECISION VHC ! HC strength
DOUBLE PRECISION PARA(10)
INTEGER SIMTYPE ! Simulation method (WLC=1,SSWLC=2,GC=3)

! Variables used for the Brownian forces

DOUBLE PRECISION FRAND(NT,3) ! Random force
DOUBLE PRECISION TRAND(NT,3) ! Random force
DOUBLE PRECISION MAGR,MAGU ! Mag of Brownian forces
INTEGER BROWN ! Logic for BD forces
INTEGER INTON ! Include polymer interactions
REAL ran1 ! Random number generator
INTEGER IDUM ! Seed for the generator
INTEGER NOW(3) ! Time now (hr,min,sec)

! Variables for the timestep switch

INTEGER SWDT

! Variable to hold time of first collisions between each bead
DOUBLE PRECISION HAS_COLLIDED(NT,NT)
DOUBLE PRECISION FPT_DIST ! l1 dist to trigger collision
INTEGER COL_TYPE ! algorithm to use for collision detection

! Variables for tracking methylation profile
INTEGER IN_RXN_RAD(NT,NT) ! is pair of sites within reaction radius? 1 = yes, 0 = no
INTEGER METH_STATUS(NT) ! methylation status of each site: 1 = methylated, 0 = unmethylated
INTEGER COULD_REACT ! number of pairs that meet criteria for reaction (i.e. within reaction radius, one bead methylated, one bead unmethylated)
INTEGER PAIRS(2,NT) ! array that holds indices of sites that could react
INTEGER RXN_HAPPEN ! reaction status: 1 = reaction, 0 = no reaction
DOUBLE PRECISION KM ! rate of methylation
DOUBLE PRECISION KD ! rate of demethylation
DOUBLE PRECISION KTOT ! total rate constant
INTEGER NUM_METHYLATED ! number of methylated sites
INTEGER NUM_SPREAD ! total number of spreading events
INTEGER NUC_SITE ! index of nucleation site
INTEGER NUM_DECAY ! total number of decay events
DOUBLE PRECISION DT_MOD ! time remaining in timestep for Gillespie algorithm

! Load the input parameters

EB=PARA(1)
EPAR=PARA(2)
EPERP=PARA(3)
GAM=PARA(4)
ETA=PARA(5)
XIR=PARA(6)
XIU=PARA(7)
LBOX=PARA(8)
LHC=PARA(9)
VHC=PARA(10)

MAGR=sqrt(XIR*2.0/DT)
MAGU=sqrt(XIU*2.0/DT)
DT0=DT
SWDT=0

! Setup the geometric parameters and initialize random forces

IB=1
DO 10 I=1,NP
DO 20 J=1,N
RS(IB,1)=R(IB,1)
RS(IB,2)=R(IB,2)
RS(IB,3)=R(IB,3)
US(IB,1)=U(IB,1)
US(IB,2)=U(IB,2)
US(IB,3)=U(IB,3)
FELAS(IB,1)=0.
FELAS(IB,2)=0.
FELAS(IB,3)=0.
FRAND(IB,1)=0.
FRAND(IB,2)=0.
FRAND(IB,3)=0.
FPONP(IB,1)=0.
FPONP(IB,2)=0.
FPONP(IB,3)=0.
TELAS(IB,1)=0.
TELAS(IB,2)=0.
TELAS(IB,3)=0.
TRAND(IB,1)=0.
TRAND(IB,2)=0.
TRAND(IB,3)=0.
TPONP(IB,1)=0.
TPONP(IB,2)=0.
TPONP(IB,3)=0.
IB=IB+1
20 CONTINUE
10 CONTINUE

RXN_HAPPEN = 1

! Begin the time integration

DO WHILE (TIME.LT.TTOT)

call CHECK_COLLISIONS(R, NT, HAS_COLLIDED, FPT_DIST, TIME, COL_TYPE, IN_RXN_RAD)

DT_MOD = DT

DO WHILE (RXN_HAPPEN.EQ.1)

COULD_REACT = 0

call CHECK_REACTIONS(R, NT, METH_STATUS, IN_RXN_RAD, COULD_REACT, FPT_DIST, PAIRS)

call TOT_RATE_CONSTANT(NT, COULD_REACT, METH_STATUS, KM, KD, KTOT, NUM_METHYLATED)

call METHYL_PROFILE(NT,METH_STATUS,KTOT,KM,KD,NUM_METHYLATED,TIME, &
RXN_HAPPEN,PAIRS,DT,DT_MOD,NUC_SITE,NUM_SPREAD,NUM_DECAY)

END DO

RXN_HAPPEN = 1

! Calculate the random forces and torques for use in this
! timestep calculation if BROWN=1

RK=1
DO WHILE (RK.LE.4)

130 CONTINUE

if (BROWN.EQ.1.AND.RK.EQ.1) then
IB=1
DO 30 I=1,NP
DO 40 J=1,N
FRAND(IB,1)=MAGR*rnorm()
FRAND(IB,2)=MAGR*rnorm()
FRAND(IB,3)=MAGR*rnorm()
TRAND(IB,1)=MAGU*rnorm()
TRAND(IB,2)=MAGU*rnorm()
TRAND(IB,3)=MAGU*rnorm()
IB=IB+1
40 CONTINUE
30 CONTINUE
endif

! Calculate the four Runge-Kutta derivatives


! Calculate the elastic forces (same as free chain)

call force_elas(FELAS,TELAS,R,U,NT,N,NP,EB,EPAR,EPERP,GAM,ETA,SIMTYPE)

! Calculate the self forces

if (INTON.EQ.1) then
call force_ponp(FPONP,R,NT,N,NP,LHC,VHC,LBOX,GAM,DT,XIR,SWDT)

! If timestep is switch, reset coords and redo step

if (SWDT.EQ.1) then
print*, "Time-step switch", DT,RK,TIME
SWDT=0
MAGR=sqrt(XIR*2.0/DT)
MAGU=sqrt(XIU*2.0/DT)
IB=1
DO 60 I=1,NP
DO 65 J=1,N
R(IB,1)=RS(IB,1)
R(IB,2)=RS(IB,2)
R(IB,3)=RS(IB,3)
U(IB,1)=US(IB,1)
U(IB,2)=US(IB,2)
U(IB,3)=US(IB,3)
IB=IB+1
65 CONTINUE
60 CONTINUE
RK=1
goto 130
endif
endif


! Calculate the change in the position vector

IB=1
DO 70 I=1,NP
DO 80 J=1,N
DRDT(IB,1,RK)=(FELAS(IB,1)+FPONP(IB,1))/XIR
DRDT(IB,2,RK)=(FELAS(IB,2)+FPONP(IB,2))/XIR
DRDT(IB,3,RK)=(FELAS(IB,3)+FPONP(IB,3))/XIR
DUDT(IB,1,RK)=(TELAS(IB,1)+TPONP(IB,1))/XIU
DUDT(IB,2,RK)=(TELAS(IB,2)+TPONP(IB,2))/XIU
DUDT(IB,3,RK)=(TELAS(IB,3)+TPONP(IB,3))/XIU

if (BROWN.EQ.0) then
DOTU=DUDT(IB,1,RK)*U(IB,1)+DUDT(IB,2,RK)*U(IB,2)+DUDT(IB,3,RK)*U(IB,3)
DUDT(IB,1,RK)=DUDT(IB,1,RK)-DOTU*U(IB,1)
DUDT(IB,2,RK)=DUDT(IB,2,RK)-DOTU*U(IB,2)
DUDT(IB,3,RK)=DUDT(IB,3,RK)-DOTU*U(IB,3)
endif
IB=IB+1
80 CONTINUE
70 CONTINUE

if (BROWN.EQ.1) then
IB=1
DO 90 I=1,NP
DO 100 J=1,N
DRDT(IB,1,RK)=DRDT(IB,1,RK)+FRAND(IB,1)/XIR
DRDT(IB,2,RK)=DRDT(IB,2,RK)+FRAND(IB,2)/XIR
DRDT(IB,3,RK)=DRDT(IB,3,RK)+FRAND(IB,3)/XIR
DUDT(IB,1,RK)=DUDT(IB,1,RK)+TRAND(IB,1)/XIU
DUDT(IB,2,RK)=DUDT(IB,2,RK)+TRAND(IB,2)/XIU
DUDT(IB,3,RK)=DUDT(IB,3,RK)+TRAND(IB,3)/XIU

DOTU=DUDT(IB,1,RK)*U(IB,1)+DUDT(IB,2,RK)*U(IB,2)+DUDT(IB,3,RK)*U(IB,3)
DUDT(IB,1,RK)=DUDT(IB,1,RK)-DOTU*U(IB,1)
DUDT(IB,2,RK)=DUDT(IB,2,RK)-DOTU*U(IB,2)
DUDT(IB,3,RK)=DUDT(IB,3,RK)-DOTU*U(IB,3)

IB=IB+1
100 CONTINUE
90 CONTINUE
endif

! If SIMTYPE=1 (WLC), calculate the constraint forces

if (SIMTYPE.EQ.1) then
call concalc(R,DRDT,NT,N,NP,XIR,GAM,DT,RK,BROWN)
endif

! Step forward using the RK algorithm

call RKstep(RS,R,US,U,DRDT,DUDT,NT,N,NP,RK,DT)

RK=RK+1

ENDDO

TIME=TIME+DT

! Swap old variables for new ones

DT=DT0
MAGR=sqrt(XIR*2.0/DT)
MAGU=sqrt(XIU*2.0/DT)

IB=1
DO 110 I=1,NP
R0(1)=nint(R(IB,1)/LBOX-0.5)*LBOX
R0(2)=nint(R(IB,2)/LBOX-0.5)*LBOX
R0(3)=nint(R(IB,3)/LBOX-0.5)*LBOX
DO 120 J=1,N
R(IB,1)=R(IB,1)-R0(1)
R(IB,2)=R(IB,2)-R0(2)
R(IB,3)=R(IB,3)-R0(3)
RS(IB,1)=R(IB,1)
RS(IB,2)=R(IB,2)
RS(IB,3)=R(IB,3)
US(IB,1)=U(IB,1)
US(IB,2)=U(IB,2)
US(IB,3)=U(IB,3)
IB=IB+1
120 CONTINUE
110 CONTINUE

ENDDO

RETURN
END

!---------------------------------------------------------------*
72 changes: 72 additions & 0 deletions BasicWLC/BDcode/RKstep.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
!---------------------------------------------------------------*

!
! This subroutine performs the RK step for user inputted timestep
! size
!
! Andrew Spakowitz
! Written 6-6-04

SUBROUTINE RKstep(RS,R,US,U,DRDT,DUDT,NT,N,NP,RK,DT)

PARAMETER (PI=3.141592654) ! Value of pi

DOUBLE PRECISION RS(NT,3) ! Saved bead positions
DOUBLE PRECISION R(NT,3) ! Temp bead positions
DOUBLE PRECISION US(NT,3) ! Unit tangent
DOUBLE PRECISION U(NT,3) ! Unit tangent
DOUBLE PRECISION DRDT(NT,3,4) ! Change rate of beads
DOUBLE PRECISION DUDT(NT,3,4) ! Change rate of beads
DOUBLE PRECISION DT ! Time step size
INTEGER N,NT ! Bead numbers
INTEGER RK ! RK number
INTEGER I,J,IB ! Index number
DOUBLE PRECISION MAGU

IB=1
DO 10 I=1,NP
DO 20 J=1,N
if(RK.EQ.1) then
R(IB,1)=RS(IB,1)+DT*DRDT(IB,1,RK)/2.
R(IB,2)=RS(IB,2)+DT*DRDT(IB,2,RK)/2.
R(IB,3)=RS(IB,3)+DT*DRDT(IB,3,RK)/2.
U(IB,1)=US(IB,1)+DT*DUDT(IB,1,RK)/2.
U(IB,2)=US(IB,2)+DT*DUDT(IB,2,RK)/2.
U(IB,3)=US(IB,3)+DT*DUDT(IB,3,RK)/2.
elseif(RK.EQ.2) then
R(IB,1)=RS(IB,1)+DT*DRDT(IB,1,RK)/2.
R(IB,2)=RS(IB,2)+DT*DRDT(IB,2,RK)/2.
R(IB,3)=RS(IB,3)+DT*DRDT(IB,3,RK)/2.
U(IB,1)=US(IB,1)+DT*DUDT(IB,1,RK)/2.
U(IB,2)=US(IB,2)+DT*DUDT(IB,2,RK)/2.
U(IB,3)=US(IB,3)+DT*DUDT(IB,3,RK)/2.
elseif(RK.EQ.3) then
R(IB,1)=RS(IB,1)+DT*DRDT(IB,1,RK)
R(IB,2)=RS(IB,2)+DT*DRDT(IB,2,RK)
R(IB,3)=RS(IB,3)+DT*DRDT(IB,3,RK)
U(IB,1)=US(IB,1)+DT*DUDT(IB,1,RK)
U(IB,2)=US(IB,2)+DT*DUDT(IB,2,RK)
U(IB,3)=US(IB,3)+DT*DUDT(IB,3,RK)
elseif(RK.EQ.4) then
R(IB,1)=RS(IB,1)+DT*(DRDT(IB,1,1)/6.+DRDT(IB,1,2)/3.+DRDT(IB,1,3)/3.+DRDT(IB,1,4)/6.)
R(IB,2)=RS(IB,2)+DT*(DRDT(IB,2,1)/6.+DRDT(IB,2,2)/3.+DRDT(IB,2,3)/3.+DRDT(IB,2,4)/6.)
R(IB,3)=RS(IB,3)+DT*(DRDT(IB,3,1)/6.+DRDT(IB,3,2)/3.+DRDT(IB,3,3)/3.+DRDT(IB,3,4)/6.)
U(IB,1)=US(IB,1)+DT*(DUDT(IB,1,1)/6.+DUDT(IB,1,2)/3.+DUDT(IB,1,3)/3.+DUDT(IB,1,4)/6.)
U(IB,2)=US(IB,2)+DT*(DUDT(IB,2,1)/6.+DUDT(IB,2,2)/3.+DUDT(IB,2,3)/3.+DUDT(IB,2,4)/6.)
U(IB,3)=US(IB,3)+DT*(DUDT(IB,3,1)/6.+DUDT(IB,3,2)/3.+DUDT(IB,3,3)/3.+DUDT(IB,3,4)/6.)
endif

MAGU=sqrt(U(IB,1)**2.+U(IB,2)**2.+U(IB,3)**2.)
U(IB,1)=U(IB,1)/MAGU
U(IB,2)=U(IB,2)/MAGU
U(IB,3)=U(IB,3)/MAGU

IB=IB+1
20 CONTINUE
10 CONTINUE

RETURN
END

!---------------------------------------------------------------*

Loading