I am facing a difficult to debug problem in this Fortran 90 code - keeping or commenting out write(*,*) statements is changing the output of the code! Here are the details and the code (cannot minimize more as that may not cause the same bug):
Os: Ubuntu 16.04
gfortran --version
GNU Fortran (Ubuntu 5.4.0-6ubuntu1~16.04.12) 5.4.0 20160609
Copyright (C) 2015 Free Software Foundation, Inc.
FORTRAN90 code name: gillsp_notworking.f90
compiling with: gfortran -ffree-line-length-none gillsp_notworking.f90
These are the two problematic lines (marked as !***** error write ****):
W1: write statement at line number 421
W2: write statement at line number 422
The code:
program code
implicit none
real*8, parameter :: time=3600.d0, t0=0.d0
real*8, parameter :: er=1.d-6
real*8, parameter :: vton=600.d0
real*8, parameter :: dl=1.d0/370.d0
real*8,parameter :: kp1be = 11.6d0/vton, km1be = 1.4d0
real*8,parameter :: kp2be = 3.4d0/vton, km2be = 0.2d0
real*8,parameter :: kp3be = 2.9d0/vton, km3be = 5.4d0
real*8,parameter :: kp1pe = 0.d0/vton, km1pe = 0.d0
real*8,parameter :: kp2pe = 0.d0/vton, km2pe = 0.d0
real*8,parameter :: kp3pe = 0.d0/vton, km3pe = 0.d0
real*8,parameter :: kpdom = 16.d-5/vton
real*8,parameter :: kpcf = 13.d0/vton, kmcf = 0.7d0
real*8,parameter :: kmbcf = 4.4, kmpcf = 0.d0
real*8,parameter :: ktd = 0.003d0
real*8,parameter :: ktdpi = 0.3d0
real*8,parameter :: ksever = 0.01d0
real*8,parameter :: kcp0=12.8d0/vton, kcm0=2.d-4
real*8,parameter :: kcp1=0.3d0/vton, kcm1=6.d-3
real*8,parameter :: alpha = 20.d0
real*8,parameter :: rho = 0.3d0
real*8,parameter :: rhocf = 0.8d0
real*8,parameter :: rhocp = 0.d0
real*8,parameter :: V = 1.d0
integer, parameter :: Nmin=100,Npos=nint(60.d0/dl)
integer, parameter :: N=nint(rho*vton*V)
integer, parameter :: Ncf=nint(rhocf*vton*V)
integer, parameter :: Ncp=nint(rhocp*vton*V/1000.d0)
integer, parameter :: nstep=300
integer, parameter :: ndmax=3000
integer, parameter :: minsevL=10
integer, parameter :: plotmovie=1
real*8, parameter :: delt=time/nstep
real*8 :: t,r1,r2,pl,ph
real*8 :: tau
real*8 :: beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8,beta_sev
real*8 :: beta_cp,beta_cm
real*8 :: beta_sum
real*8 :: alpha0
real*8 :: rsev1,rsev2
real*8 :: pifactr
integer :: m, m1(-Nmin:Npos), m10(-Nmin:Npos)
integer :: mcf, domL(1:ndmax),domR(1:ndmax), ndom,cofcount,domLen
integer :: cap
integer :: domL0(1:ndmax),domR0(1:ndmax),drcount,domrem(1:ndmax)
integer :: pesever(1:10),besever(1:10),besvcount,pesvcount, sevindex, dumind
integer :: i,j,k,nk
integer :: l1(-Nmin:Npos)
integer :: l1bind, l1pind, l1bind0, l1pind0
integer :: gatp1, gadppi1, gadp1, glcf1, gzero
integer :: l1sum
integer :: k1,k2, dgcount, dscount, nmerge, domflag
integer :: nfile
integer :: dumdl
integer :: Bstate
integer :: nPicof
integer :: seed(12)
character*80::fname, Rid
write(*,*) 'Array size=', -Nmin, Npos
!------------------------------------------------------
seed= 9878771
CALL RANDOM_SEED (PUT=seed)
!------------------------------------------------------
t=0
nk=0
nfile=0
!----------- initialisation -----------
l1=0
m1=0
l1bind=0
l1pind=-1
m=N
cap=Ncp
Bstate = 1
pifactr = 1.d0
ndom=0
cofcount=0
domL=0
domR=0
nPicof=0
mcf=Ncf-sum(domR-domL)-ndom
!********************************************************************
!********************************************************************
!------------------- time loop ------------------------
tloop: do while (t.le.time)
!----------------- counting G-atp, G-adp, G-cofilin in filament --------------
gzero=count(m1<er)
gatp1=count(m1<1+er) - gzero
gadppi1=count(m1<2+er) - (gatp1+gzero)
gadp1=count(m1<3+er) - (gatp1+gzero+gadppi1)
glcf1=count(m1>3+er)
!----------------- current length ---------------------
l1sum = sum(l1)
!------------------ propensity sum ---------------------
!---- store monomer states, domain info etc -----
m10 = m1
l1bind0=l1bind
l1pind0=l1pind
domL0=domL
domR0=domR
!------------------- capping reaction ------------------
if (m10(l1bind-1).eq.4) then
if(Bstate.eq.1)then
beta_cp = kcp1*cap/V
beta_cm = 0.d0
elseif (Bstate.eq.2) then
beta_cp = 0.d0
beta_cm = kcm1
endif
else
if(Bstate.eq.1)then
beta_cp = kcp0*cap/V
beta_cm = 0.d0
elseif (Bstate.eq.2) then
beta_cp = 0.d0
beta_cm = kcm0
endif
endif
!-------------------- barbed end -----------------------
if(Bstate.eq.1)then ! Barbed-end state if
if(m10(l1bind-1).eq.1)then
beta1 = km1be
beta2 = kp1be*(m/V)
elseif (m10(l1bind-1).eq.2) then
beta1 = km2be
beta2 = kp2be*(m/V)
elseif (m10(l1bind-1).eq.3) then
beta1 = km3be
beta2 = kp3be*(m/V)
elseif (m10(l1bind-1).eq.4) then
beta1 = kmbcf
beta2 = 0.d0
endif
elseif (Bstate.eq.2) then
beta1 = 0.d0
beta2 = 0.d0
endif
!----- the first monomer incorporation------!
if(t.lt.30)then
if(m10(l1bind-1).eq.0)then
beta2 = kp1be*(m/V)
endif
endif
!-------------------- pointed end -----------------------
if(m10(l1pind+1).eq.1) then
beta3 = km1pe
beta4=kp1pe*(m/V)
elseif (m10(l1pind+1).eq.2) then
beta3 = km2pe
beta4=kp2pe*(m/V)
elseif (m10(l1pind+1).eq.3) then
beta3 = km3pe
beta4 = kp3pe*(m/V)
elseif (m10(l1pind+1).eq.4) then
beta3 = kmpcf
beta4 = 0.d0
endif
!------------------- bulk reactions ---------------------
!------------------ cofilin reactions ------------------
if(ndom.gt.0)then
nmerge=0
do k1=1,ndom
do k2=1,ndom
if(domL(k1)-1.eq.domR(k2)) nmerge=nmerge+1
enddo
enddo
endif
!---------------- Pi-cofilin interface ------------------
nPicof=0
do k1=l1pind0,l1bind0-1
if(abs(m10(k1)-m10(k1+1)).eq.2.and.m10(k1)+m10(k1+1).eq.6)then
nPicof=nPicof+1
endif
enddo
!------------------- bulk reactions ---------------------
beta5 = ktdpi*gatp1 + ktd*(gadppi1-nPicof) + alpha*ktd*nPicof
beta6 = kpdom*(mcf/V)*max(0.d0,(gadp1-2.d0*ndom-2.d0*nmerge))
!----- count domain growth probability -----
if(ndom.gt.0)then
dgcount=0
dscount=0
do k=1,ndom
dumdl=domR(k)-domL(k)-1
if(m10(domL(k)).eq.3) dgcount=dgcount+1
if(m10(domR(k)).eq.3) dgcount=dgcount+1
if(m10(domL(k)).le.3.and.dumdl.gt.2) dscount=dscount+1
if(m10(domR(k)).le.3.and.dumdl.gt.2) dscount=dscount+1
enddo
else
dgcount=0
dscount=0
endif
beta7 = kpcf*(mcf/V)*dgcount
beta8 = kmcf*dscount
!----- severing probability -----
if(ndom.gt.0)then
besvcount=0
pesvcount=0
besever=-1
pesever=-1
do k=1,ndom
domLen = domR(k)-domL(k)-1
if (domLen.ge.minsevL) then
if(m10(domL(k)).lt.4.and.m10(domL(k)).gt.0) then
pesvcount=pesvcount+1
pesever(pesvcount)=domL(k)
endif
if(m10(domR(k)).lt.4.and.m10(domR(k)).gt.0) then
besvcount=besvcount+1
besever(besvcount)=domR(k)
endif
else !------ domain length condition
pesvcount=0
besvcount=0
endif !------ domain length condition
enddo
else
pesvcount=0
besvcount=0
endif
beta_sev = ksever*(besvcount+pesvcount)
beta_sum = beta1+beta2+beta3+beta4+beta5+beta6+beta7 &
+beta8+beta_sev+beta_cp+beta_cm
alpha0 = beta_sum
!----------------- next reaction time: tau --------------
CALL RANDOM_NUMBER(r1)
if(r1.lt.1.d-8)then
CALL RANDOM_NUMBER(r1)
endif
tau = (1.d0/alpha0)*log(1.d0/r1)
!----------------------- all reactions ------------------------
CALL RANDOM_NUMBER(r2)
!---------------- growth and decay ------------------
!---------------- filament 1 ----------------
!---------------- barbed end ----------------
pl=0.d0
ph=beta1/alpha0
if (r2.ge.pl.and.r2.lt.ph.and.l1sum.gt.er) then
if (m1(l1bind-1).eq.4) then
! mcf = mcf+1 Pool conc. remains same
do k=1,ndom
if(domR(k).eq.l1bind) domR(k)=domR(k)-1
enddo
endif
l1bind = l1bind - 1
l1(l1bind) = 0
m1(l1bind) = 0
! m = m+1 Pool conc. remains same
Rid="BE-"
endif
pl=beta1/alpha0
ph=(beta1+beta2)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
m1(l1bind) = 1
l1(l1bind) = 1
l1bind = l1bind + 1
! m = m-1 Pool conc. remains same
Rid="BE+"
endif
!---------------- pointed end ---------------
pl=(beta1+beta2)/alpha0
ph=(beta1+beta2+beta3)/alpha0
if (r2.ge.pl.and.r2.lt.ph.and.l1sum.gt.er) then
if (m1(l1pind+1).eq.4) then
! mcf = mcf+1 Pool conc. remains same
do k=1,ndom
if(domL(k).eq.l1pind) domL(k)=domL(k)+1
enddo
endif
l1pind = l1pind + 1
l1(l1pind) = 0
m1(l1pind) = 0
! m = m+1 Pool conc. remains same
Rid="PE-"
endif
pl=(beta1+beta2+beta3)/alpha0
ph=(beta1+beta2+beta3+beta4)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
m1(l1pind) = 1
l1(l1pind) = 1
l1pind = l1pind - 1
! m = m-1 Pool conc. remains same
Rid="PE+"
endif
!---------- severing of filament -------------
pl=(beta1+beta2+beta3+beta4)/alpha0
ph=(beta1+beta2+beta3+beta4+beta_sev)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
CALL RANDOM_NUMBER(rsev1)
CALL RANDOM_NUMBER(rsev2)
if(rsev1.le.0.2d0.and.besvcount.gt.0)then
dumind = 1 + nint(rsev2*(besvcount-1))
sevindex = besever(dumind)
else
if (pesvcount.gt.0) then
dumind = 1 + nint(rsev2*(pesvcount-1))
sevindex = pesever(dumind)
else
dumind = 1 + nint(rsev2*(besvcount-1))
sevindex = besever(dumind)
endif
endif
l1bind=sevindex ! keep olp pointed end
Bstate=1 ! set new barbed end free
! m=m+(l1bind0-sevindex) Pool conc. remains same
m1(sevindex:l1bind0)=0
l1(sevindex:l1bind0)=0
drcount=0
domrem=0
do j=sevindex,l1bind0 ! j-loop
! if(m10(j).eq.4) then Pool conc. remains same
! mcf=mcf+1
! endif
do k=1,ndom
if(domL(k).eq.j)then
drcount=drcount+1
domrem(drcount)=j
domL(k)=0
domR(k)=0
endif
enddo
enddo ! j-loop
Rid="SEVER"
! write(*,*)'sever', m10(sevindex-1) ,sevindex, l1bind0, besvcount, pesvcount !***** error write ****
! write(*,*)'sever', m10(sevindex-1) !***** error write ****
endif ! filament severing if
!---------------- capping ------------------
pl=ph
ph=(beta1+beta2+beta3+beta4+beta_sev+beta_cp)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
! cap=cap-1 Pool conc. remains same
Bstate = 2
Rid="CAPPING +"
endif
!---------------- de-capping ------------------
pl=ph
ph=(beta1+beta2+beta3+beta4+beta_sev+beta_cp+beta_cm)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
! cap=cap+1 Pool conc. remains same
Bstate = 1
Rid="CAPPING -"
endif
!---------- hydrolysis in filament -----------
pl=(beta1+beta2+beta3+beta4+beta_sev+beta_cp+beta_cm)/alpha0
do j=l1pind0,l1bind0
if(m10(j).eq.1)then
ph=pl+(ktdpi/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 2
Rid="HYD atp"
endif
pl=ph
elseif(m10(j).eq.2)then
if(m10(j-1).eq.4.or.m10(j+1).eq.4)then
pifactr=alpha
else
pifactr=1.d0
endif
ph=pl+(pifactr*ktd/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 3
Rid="Pi Release"
endif
pl=ph
elseif(m10(j).eq.3)then
domflag=0
if(ndom.gt.0)then ! ndom-if
do k=1,ndom
if(j.eq.domL(k)) then
domflag=1
endif
if(j.eq.domR(k)) then
domflag=1
endif
enddo
endif ! ndom-if
if (domflag.lt.1) then
ph=pl+((kpdom*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 4
! mcf=mcf-1 Pool conc. remains same
ndom=ndom+1
domL(ndom)=j-1
domR(ndom)=j+1
Rid="DOM"
endif
pl=ph
endif
endif ! hydrolysis if
enddo ! j monomer loop
if(ndom.gt.0)then ! ndom if
do k=1,ndom
dumdl=domR(k)-domL(k)-1
!--------------------- domain growth --------------------
if (m10(domL(k)).eq.3) then ! Left boundary ADP-if
ph=pl+((kpcf*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domL(k))=4
! mcf=mcf-1 Pool conc. remains same
domL(k)=domL(k)-1
Rid="DOML+"
endif
pl=ph
endif
if (m10(domR(k)).eq.3) then ! Right boundary ADP-if
ph=pl+((kpcf*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domR(k))=4
! mcf=mcf-1 Pool conc. remains same
domR(k)=domR(k)+1
Rid="DOMR+"
endif
pl=ph
endif
! --------------------- domain shrinkage -------------------
if (m10(domL(k)).le.3.and.dumdl.gt.2) then
ph=pl+(kmcf/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domL(k)+1)=3
! mcf=mcf+1 Pool conc. remains same
domL(k)=domL(k)+1
Rid="DOML-"
endif
pl=ph
endif
if (m10(domR(k)).le.3.and.dumdl.gt.2) then
ph=pl+(kmcf/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domR(k)-1)=3
! mcf=mcf+1 Pool conc. remains same
domR(k)=domR(k)-1
Rid="DOMR-"
endif
pl=ph
endif
enddo ! k ndom loop
endif ! ndom if
!----------------- count bound cofilin --------------------
cofcount=0
do j=l1pind,l1bind
if (m1(j).eq.4) then
cofcount=cofcount+1
endif
enddo
!------------------- ERROR CHECKS ----------------------
if (tau.lt.0.d0) exit tloop
if (m1(l1bind).ne.0.or.m1(l1pind).ne.0) exit tloop
if (m.ne.N) exit tloop
if (mcf.ne.Ncf) exit tloop
if (cap.ne.Ncp) exit tloop
t = t + tau
enddo tloop !***** time loop *****
!------------------ time loop ended -----------------------
write(*,*) 'domains created', ndom
do k=1,ndom
write(*,*) k, 'size', domR(k)-domL(k)-1, 'L-R', domL(k), domR(k)
enddo
do k=l1pind,l1bind
write(101,*) k, m1(k)
enddo
if(t.ge.time) write(*,*) 'CODE ENDED WITHOUT ERROR'
write(*,*) 'CODE STOPPED ','t=',t,'id=',Rid
if (tau.lt.0.d0) write(*,*) 'ERROR: tau negative', tau, alpha0, beta1+beta2+beta3+beta4+beta5, beta6, beta7
if (m1(l1bind).ne.0.or.m1(l1pind).ne.0) write(*,*) 'ERROR: l1bind/l1pind not zero', m1(l1bind), m1(l1pind)
if (m.ne.N) write(*,*) 'ERROR: constant conc. fails', m, N
if (mcf.ne.Ncf) write(*,*) 'ERROR: constant conc. fails', mcf, Ncf
if (cap.ne.Ncp) write(*,*) 'ERROR: constant conc. fails', cap, Ncp
write(*,*) 'pe',l1pind,'be',l1bind, 'domains', ndom
stop
end program code
O/P when W1 and W2 are commented out:
Array size= -100 22200
domains created 30
1 size 40 L-R 0 41
2 size 0 L-R -1 0
3 size -1 L-R 0 0
4 size -1 L-R 0 0
5 size -1 L-R 0 0
6 size -1 L-R 0 0
7 size 113 L-R 40 154
8 size -1 L-R 0 0
9 size -1 L-R 0 0
10 size -1 L-R 0 0
11 size -1 L-R 0 0
12 size 57 L-R 153 211
13 size -40 L-R 250 211
14 size -1 L-R 0 0
15 size -1 L-R 0 0
16 size -11 L-R 378 368
17 size -1 L-R 0 0
18 size -1 L-R 0 0
19 size -1 L-R 0 0
20 size -1 L-R 0 0
21 size -1 L-R 0 0
22 size -16 L-R 499 484
23 size -1 L-R 0 0
24 size -8 L-R 545 538
25 size -5 L-R 581 577
26 size -11 L-R 516 506
27 size -41 L-R 617 577
28 size -1 L-R 0 0
29 size -1 L-R 0 0
30 size -1 L-R 703 703
CODE ENDED WITHOUT ERROR
CODE STOPPED t= 3600.9739030733977 id=BE+
pe -1 be 239 domains 30
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL
O/P when W1 is NOT commented out:
Array size= -100 22200
sever 4 144 550 2 1
sever 2 108 239 2 1
sever 2 207 475 1 1
sever 2 256 407 2 2
sever 2 224 579 2 2
sever 4 226 297 1 1
sever 2 337 562 2 2
sever 3 324 476 1 1
sever 4 376 611 1 1
sever 4 369 436 1 0
sever 4 615 943 1 1
sever 3 582 776 1 2
sever 4 409 734 1 1
sever 2 364 403 1 1
sever 4 512 833 3 2
sever 3 463 568 2 2
sever 4 450 493 1 1
sever 2 430 446 0 1
sever 3 865 1299 4 4
sever 4 812 1041 3 3
sever 4 547 813 1 1
sever 2 545 851 2 1
sever 2 573 765 1 1
domains created 31
1 size 122 L-R 0 123
2 size -1 L-R 0 0
3 size 0 L-R -1 0
4 size -1 L-R 0 0
5 size 75 L-R 122 198
6 size -1 L-R 0 0
7 size -1 L-R 0 0
8 size 39 L-R 197 237
9 size 48 L-R 236 285
10 size -1 L-R 0 0
11 size 3 L-R 316 320
12 size -1 L-R 0 0
13 size 0 L-R 356 357
14 size -1 L-R 0 0
15 size 36 L-R 319 356
16 size -1 L-R 0 0
17 size -1 L-R 0 0
18 size 139 L-R 357 497
19 size -1 L-R 0 0
20 size -1 L-R 0 0
21 size -1 L-R 0 0
22 size -1 L-R 0 0
23 size -1 L-R 0 0
24 size -1 L-R 0 0
25 size -1 L-R 0 0
26 size -1 L-R 0 0
27 size -1 L-R 0 0
28 size 47 L-R 517 565
29 size 21 L-R 496 518
30 size -1 L-R 0 0
31 size -1 L-R 0 0
CODE ENDED WITHOUT ERROR
CODE STOPPED t= 3600.0320511744094 id=HYD atp
pe -1 be 688 domains 31
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL
None of the results are unphysical and the different results can be reproduced with a different compiler on a mac. I have also compiled with gfortran -Wall -Wno-tabs -g -fcheck=all -fbacktrace gillsp_notworking.f90 but did not find any serious issue. Apart from the different results, the code can sometime freeze as well in some combination of commenting out the two write statements and that seems to depend on what combination I ran before (crazy!) - if this helps in the diagnosis. Any help or hints are welcome. Thank you in advance!
------ EDIT 1 --------
compiling with optimization flag supresses different outputs but freezes in some cases (e.g., when W1,W2 are not commented out):
gfortran -O gillsp_notworking.f90
does not solve the problem.
------ EDIT 2 ---------
Following suggestion from comments, I have now compiled with debugging options and bound check options:
gfortran -Wall -Wno-tabs -Wextra -fbounds-check gillsp_notworking.f90
but this just gives me some warning about unused parameters.
Following the suggestion of using random_number() I have modified the code in the current version and the subsequent different behaviours too. It freezes if I keep both the write statements. Also if I compile using -Og flag (i.e., gfortran -Og gillsp_notworking_seed.f90) when both W1,W2 are commented out then the results change!
Using GFortran 11.3 (the default version included in Ubuntu 22.04), compiling with -Og -Wall and then fixing the may be used uninitialized warnings by setting in the initialization block
beta1 = 0
beta2 = 0
beta3 = 0
beta4 = 0
beta5 = 0
beta_cp = 0
beta_cm = 0
nmerge = 0
I get the same output regardless if those two write statements are commented out or not, and regardless of optimization level (tried -O0, -O2, -Og, -O3 and -O3 -ffast-math). Also neither -fcheck=all and -fsanitize=undefined nor -fsanitize=address finds any problem (which of course isn't a guarantee that problems aren't there, but gives a bit of indication at least).