!======================================================================= ! NVT Monte Carlo Simulation — Main Program ! All parameters are read from 'input.macro' at runtime. ! DO NOT hardcode any parameters here. !======================================================================= implicit none !--- Integer variables --- integer n,updateInter,nStep,i,j,atomNo,nPrint,printCount,w1,v1,max_attempts integer outputUnit,bend_samples,flexibility_output_freq,equilFrac,xyzFreq !--- Double precision variables --- double precision temp,rCut,beta,sigma,rMin,maxDisplacement double precision x(2000,10),y(2000,10),z(2000,10),v,vRecalc double precision ex(2000,10),ey(2000,10),ez(2000,10) double precision vEnd,vNew,vOld,a,b,c double precision delBV,delV,prob,acceptCount,stepCount double precision ratio,box,boxInverse,bond_length_tail_tail,volume double precision x1Old,y1Old,z1Old,x2Old,y2Old,z2Old,x3Old,y3Old,z3Old double precision x4Old,y4Old,z4Old,x5Old,y5Old,z5Old,x6Old,y6Old,z6Old double precision x7Old,y7Old,z7Old,x8Old,y8Old,z8Old,x9Old,y9Old,z9Old double precision x10Old,y10Old,z10Old double precision x1New,y1New,z1New,x2New,y2New,z2New,x3New,y3New,z3New double precision x4New,y4New,z4New,x5New,y5New,z5New,x6New,y6New,z6New double precision x7New,y7New,z7New,x8New,y8New,z8New,x9New,y9New,z9New double precision x10New,y10New,z10New double precision acceptanceRatio,acceptance,rot_acceptanceRatio,rot_acceptance,rot_ratio double precision dphimax,dcosmax,pi double precision rot_acceptCount,apply_pbc double precision avg_end_to_end,avg_rg,sum_end_to_end,sum_rg,sum_end_to_end_sq,sum_rg_sq double precision avg_bond_angle,sum_bond_angles(9),persistence_length double precision wF,wCut,attraction,frac,w11_param,kappa_param !--- Logical --- logical overlap !--- Filenames --- character(len=256) inputFile,outputDat,outputXyz,outputFlex character(len=256) macroFile,keyword,line,strval !--- Common blocks --- common /position1/x,y,z,ex,ey,ez,n,atomNo common /flexibility_stats/ sum_end_to_end,sum_rg,sum_end_to_end_sq,sum_rg_sq,sum_bond_angles,bend_samples common /potparams/ wF,wCut,attraction,frac,w11_param,kappa_param !======================================================================= ! READ ALL PARAMETERS FROM input.macro !======================================================================= macroFile = 'input.macro' outputUnit = 101 ! --- Defaults (fallback if any keyword missing) --- n=100; atomNo=10; nStep=50000; updateInter=50; xyzFreq=1000 temp=1.0D0; box=23.0D0; sigma=1.0D0; rMin=0.7D0 wF=0.5D0; wCut=2.5D0; attraction=0.5D0; frac=0.2D0 w11_param=8.0D0; kappa_param=1.0D0/7.0D0 maxDisplacement=0.5D0; dPhiMax=0.1D0; dCosMax=0.1D0; max_attempts=40 flexibility_output_freq=1000; equilFrac=10 inputFile = 'input.dat' outputDat = 'output.dat' outputXyz = 'output.xyz' outputFlex = 'flexibility_analysis.dat' open(unit=50, file=macroFile, status='old', err=900) do read(50,'(A)',end=10) line line = adjustl(line) if (len_trim(line)==0) cycle if (line(1:1)=='!') cycle read(line,*,err=5) keyword, strval select case(trim(keyword)) case('n') ; read(strval,*) n case('atomNo') ; read(strval,*) atomNo case('nStep') ; read(strval,*) nStep case('updateInter') ; read(strval,*) updateInter case('xyzFreq') ; read(strval,*) xyzFreq case('temp') ; read(strval,*) temp case('box') ; read(strval,*) box case('sigma') ; read(strval,*) sigma case('rMin') ; read(strval,*) rMin case('wF') ; read(strval,*) wF case('wCut') ; read(strval,*) wCut case('attraction') ; read(strval,*) attraction case('frac') ; read(strval,*) frac case('w11') ; read(strval,*) w11_param case('kappa') ; read(strval,*) kappa_param case('maxDisp') ; read(strval,*) maxDisplacement case('dPhiMax') ; read(strval,*) dPhiMax case('dCosMax') ; read(strval,*) dCosMax case('maxAttempts') ; read(strval,*) max_attempts case('flexFreq') ; read(strval,*) flexibility_output_freq case('equilFrac') ; read(strval,*) equilFrac case('inputFile') ; read(line,*,err=5) keyword, inputFile case('outputDat') ; read(line,*,err=5) keyword, outputDat case('outputXyz') ; read(line,*,err=5) keyword, outputXyz case('outputFlex') ; read(line,*,err=5) keyword, outputFlex case default write(*,*) "WARNING: Unknown macro keyword ignored: ", trim(keyword) end select 5 continue enddo 10 close(50) goto 20 900 write(*,*) "ERROR: Cannot open macro file: ", trim(macroFile); stop 20 continue !======================================================================= ! DERIVED PARAMETERS !======================================================================= pi = 4.0D0*atan(1.0D0) beta = 1.0D0/temp boxInverse = 1.0D0/box volume = box**3 rMin = rMin*sigma wF = wF*sigma wCut = wCut*sigma bond_length_tail_tail = 1.0D0*sigma maxDisplacement = maxDisplacement*sigma nPrint = nStep/100 printCount = 1 stepCount = 0.0D0 v = 0.0D0 acceptCount = 0.0D0 rot_acceptCount = 0.0D0 acceptance = 0.0D0 rot_acceptance = 0.0D0 rCut = 2.5D0*sigma if (rCut > 0.5D0*box) rCut = 0.5D0*box bend_samples = 0 sum_end_to_end=0.0D0; sum_rg=0.0D0 sum_end_to_end_sq=0.0D0; sum_rg_sq=0.0D0 do i=1,9; sum_bond_angles(i)=0.0D0; enddo !======================================================================= ! ECHO ALL PARAMETERS !======================================================================= write(*,*) "======================================================" write(*,*) " NVT MC SIMULATION — parameters read from input.macro" write(*,*) "======================================================" write(*,*) "Input file : ", trim(inputFile) write(*,*) "Output dat file : ", trim(outputDat) write(*,*) "Output xyz file : ", trim(outputXyz) write(*,*) "Output flex file : ", trim(outputFlex) write(*,*) "n (molecules) :", n write(*,*) "atomNo (beads/mol) :", atomNo write(*,*) "nStep :", nStep write(*,*) "XYZ write freq :", xyzFreq write(*,*) "updateInter :", updateInter write(*,*) "temp :", temp write(*,*) "box :", box write(*,*) "sigma :", sigma write(*,*) "rMin :", rMin write(*,*) "wF :", wF write(*,*) "wCut :", wCut write(*,*) "attraction :", attraction write(*,*) "frac :", frac write(*,*) "w11 :", w11_param write(*,*) "kappa :", kappa_param write(*,*) "maxDisp :", maxDisplacement write(*,*) "dPhiMax :", dPhiMax write(*,*) "dCosMax :", dCosMax write(*,*) "max_attempts :", max_attempts write(*,*) "flexFreq :", flexibility_output_freq write(*,*) "equilFrac :", equilFrac write(*,*) "======================================================" !======================================================================= ! READ INITIAL CONFIGURATION !======================================================================= open(unit=10, file=trim(inputFile), status='old', err=901) do i=1,n do j=1,atomNo read(10,*) x(i,j),y(i,j),z(i,j) x(i,j)=apply_pbc(x(i,j),box) y(i,j)=apply_pbc(y(i,j),box) z(i,j)=apply_pbc(z(i,j),box) enddo enddo close(10) goto 21 901 write(*,*) "ERROR: Cannot open input file: ", trim(inputFile); stop 21 continue !======================================================================= ! INITIAL ENERGY & VALIDATION !======================================================================= if (rCut .gt. 0.5D0*box) then write(*,*) "ERROR: rCut too large for box size"; stop endif call potentialCal(box,sigma,rMin,temp,overlap,v) write(*,*) "Initial Potential =", v if (overlap) then write(*,*) "ERROR: Overlap in initial configuration"; stop endif ! Create fresh (empty) output files open(unit=outputUnit, file=trim(outputDat), status='replace'); close(outputUnit) open(unit=outputUnit+1, file=trim(outputXyz), status='replace'); close(outputUnit+1) !======================================================================= ! MAIN MC LOOP !======================================================================= do i=1,nStep if (i .eq. (nPrint*printCount)) then write(*,*) ((i*100)/nStep),"%" printCount=printCount+1 endif do j=1,n ! Store old positions x1Old=x(j,1); y1Old=y(j,1); z1Old=z(j,1) x2Old=x(j,2); y2Old=y(j,2); z2Old=z(j,2) x3Old=x(j,3); y3Old=y(j,3); z3Old=z(j,3) x4Old=x(j,4); y4Old=y(j,4); z4Old=z(j,4) x5Old=x(j,5); y5Old=y(j,5); z5Old=z(j,5) x6Old=x(j,6); y6Old=y(j,6); z6Old=z(j,6) x7Old=x(j,7); y7Old=y(j,7); z7Old=z(j,7) x8Old=x(j,8); y8Old=y(j,8); z8Old=z(j,8) x9Old=x(j,9); y9Old=y(j,9); z9Old=z(j,9) x10Old=x(j,10); y10Old=y(j,10); z10Old=z(j,10) call energyCal(x1Old,y1Old,z1Old,x2Old,y2Old,z2Old,x3Old,y3Old,z3Old, & x4Old,y4Old,z4Old,x5Old,y5Old,z5Old,x6Old,y6Old,z6Old, & x7Old,y7Old,z7Old,x8Old,y8Old,z8Old,x9Old,y9Old,z9Old, & x10Old,y10Old,z10Old,j,box,sigma,temp,vOld) !--- TRANSLATION --- call random_number(a); call random_number(b); call random_number(c) x1New=apply_pbc(x1Old+((2.0D0*a-1)*maxDisplacement),box) y1New=apply_pbc(y1Old+((2.0D0*b-1)*maxDisplacement),box) z1New=apply_pbc(z1Old+((2.0D0*c-1)*maxDisplacement),box) x2New=apply_pbc(x2Old+((2.0D0*a-1)*maxDisplacement),box) y2New=apply_pbc(y2Old+((2.0D0*b-1)*maxDisplacement),box) z2New=apply_pbc(z2Old+((2.0D0*c-1)*maxDisplacement),box) x3New=apply_pbc(x3Old+((2.0D0*a-1)*maxDisplacement),box) y3New=apply_pbc(y3Old+((2.0D0*b-1)*maxDisplacement),box) z3New=apply_pbc(z3Old+((2.0D0*c-1)*maxDisplacement),box) x4New=apply_pbc(x4Old+((2.0D0*a-1)*maxDisplacement),box) y4New=apply_pbc(y4Old+((2.0D0*b-1)*maxDisplacement),box) z4New=apply_pbc(z4Old+((2.0D0*c-1)*maxDisplacement),box) x5New=apply_pbc(x5Old+((2.0D0*a-1)*maxDisplacement),box) y5New=apply_pbc(y5Old+((2.0D0*b-1)*maxDisplacement),box) z5New=apply_pbc(z5Old+((2.0D0*c-1)*maxDisplacement),box) x6New=apply_pbc(x6Old+((2.0D0*a-1)*maxDisplacement),box) y6New=apply_pbc(y6Old+((2.0D0*b-1)*maxDisplacement),box) z6New=apply_pbc(z6Old+((2.0D0*c-1)*maxDisplacement),box) x7New=apply_pbc(x7Old+((2.0D0*a-1)*maxDisplacement),box) y7New=apply_pbc(y7Old+((2.0D0*b-1)*maxDisplacement),box) z7New=apply_pbc(z7Old+((2.0D0*c-1)*maxDisplacement),box) x8New=apply_pbc(x8Old+((2.0D0*a-1)*maxDisplacement),box) y8New=apply_pbc(y8Old+((2.0D0*b-1)*maxDisplacement),box) z8New=apply_pbc(z8Old+((2.0D0*c-1)*maxDisplacement),box) x9New=apply_pbc(x9Old+((2.0D0*a-1)*maxDisplacement),box) y9New=apply_pbc(y9Old+((2.0D0*b-1)*maxDisplacement),box) z9New=apply_pbc(z9Old+((2.0D0*c-1)*maxDisplacement),box) x10New=apply_pbc(x10Old+((2.0D0*a-1)*maxDisplacement),box) y10New=apply_pbc(y10Old+((2.0D0*b-1)*maxDisplacement),box) z10New=apply_pbc(z10Old+((2.0D0*c-1)*maxDisplacement),box) call random_number(prob) call test(x1New,y1New,z1New,x2New,y2New,z2New,x3New,y3New,z3New, & x4New,y4New,z4New,x5New,y5New,z5New,x6New,y6New,z6New, & x7New,y7New,z7New,x8New,y8New,z8New,x9New,y9New,z9New, & x10New,y10New,z10New,j,box,rMin,overlap) if (.not. overlap) then call energyCal(x1New,y1New,z1New,x2New,y2New,z2New,x3New,y3New,z3New, & x4New,y4New,z4New,x5New,y5New,z5New,x6New,y6New,z6New, & x7New,y7New,z7New,x8New,y8New,z8New,x9New,y9New,z9New, & x10New,y10New,z10New,j,box,sigma,temp,vNew) delV=vNew-vOld; delBV=beta*delV if (delBV .lt. 75.0D0) then if (delBV.le.0.0D0 .or. exp(-delBV).gt.prob) then x(j,1)=x1New; y(j,1)=y1New; z(j,1)=z1New x(j,2)=x2New; y(j,2)=y2New; z(j,2)=z2New x(j,3)=x3New; y(j,3)=y3New; z(j,3)=z3New x(j,4)=x4New; y(j,4)=y4New; z(j,4)=z4New x(j,5)=x5New; y(j,5)=y5New; z(j,5)=z5New x(j,6)=x6New; y(j,6)=y6New; z(j,6)=z6New x(j,7)=x7New; y(j,7)=y7New; z(j,7)=z7New x(j,8)=x8New; y(j,8)=y8New; z(j,8)=z8New x(j,9)=x9New; y(j,9)=y9New; z(j,9)=z9New x(j,10)=x10New; y(j,10)=y10New; z(j,10)=z10New v=v+delV; acceptCount=acceptCount+1.0D0; acceptance=acceptance+1.0D0 endif endif endif stepCount=stepCount+1.0D0 ! Save accepted translation state x1New=x(j,1); y1New=y(j,1); z1New=z(j,1) x2New=x(j,2); y2New=y(j,2); z2New=z(j,2) x3New=x(j,3); y3New=y(j,3); z3New=z(j,3) x4New=x(j,4); y4New=y(j,4); z4New=z(j,4) x5New=x(j,5); y5New=y(j,5); z5New=z(j,5) x6New=x(j,6); y6New=y(j,6); z6New=z(j,6) x7New=x(j,7); y7New=y(j,7); z7New=z(j,7) x8New=x(j,8); y8New=y(j,8); z8New=z(j,8) x9New=x(j,9); y9New=y(j,9); z9New=z(j,9) x10New=x(j,10); y10New=y(j,10); z10New=z(j,10) !--- ROTATION --- call energyCal(x1New,y1New,z1New,x2New,y2New,z2New,x3New,y3New,z3New, & x4New,y4New,z4New,x5New,y5New,z5New,x6New,y6New,z6New, & x7New,y7New,z7New,x8New,y8New,z8New,x9New,y9New,z9New, & x10New,y10New,z10New,j,box,sigma,temp,vOld) call rotate_flexible_molecule(j,box,rMin,bond_length_tail_tail,max_attempts) call random_number(prob) x1Old=x(j,1); y1Old=y(j,1); z1Old=z(j,1) x2Old=x(j,2); y2Old=y(j,2); z2Old=z(j,2) x3Old=x(j,3); y3Old=y(j,3); z3Old=z(j,3) x4Old=x(j,4); y4Old=y(j,4); z4Old=z(j,4) x5Old=x(j,5); y5Old=y(j,5); z5Old=z(j,5) x6Old=x(j,6); y6Old=y(j,6); z6Old=z(j,6) x7Old=x(j,7); y7Old=y(j,7); z7Old=z(j,7) x8Old=x(j,8); y8Old=y(j,8); z8Old=z(j,8) x9Old=x(j,9); y9Old=y(j,9); z9Old=z(j,9) x10Old=x(j,10); y10Old=y(j,10); z10Old=z(j,10) call energyCal(x1Old,y1Old,z1Old,x2Old,y2Old,z2Old,x3Old,y3Old,z3Old, & x4Old,y4Old,z4Old,x5Old,y5Old,z5Old,x6Old,y6Old,z6Old, & x7Old,y7Old,z7Old,x8Old,y8Old,z8Old,x9Old,y9Old,z9Old, & x10Old,y10Old,z10Old,j,box,sigma,temp,vNew) delV=vNew-vOld; delBV=beta*delV if (delBV .lt. 75.0D0) then if (delBV.le.0.0D0 .or. exp(-delBV).gt.prob) then v=v+delV; rot_acceptCount=rot_acceptCount+1.0D0; rot_acceptance=rot_acceptance+1.0D0 else x(j,1)=x1New; y(j,1)=y1New; z(j,1)=z1New x(j,2)=x2New; y(j,2)=y2New; z(j,2)=z2New x(j,3)=x3New; y(j,3)=y3New; z(j,3)=z3New x(j,4)=x4New; y(j,4)=y4New; z(j,4)=z4New x(j,5)=x5New; y(j,5)=y5New; z(j,5)=z5New x(j,6)=x6New; y(j,6)=y6New; z(j,6)=z6New x(j,7)=x7New; y(j,7)=y7New; z(j,7)=z7New x(j,8)=x8New; y(j,8)=y8New; z(j,8)=z8New x(j,9)=x9New; y(j,9)=y9New; z(j,9)=z9New x(j,10)=x10New; y(j,10)=y10New; z(j,10)=z10New endif else x(j,1)=x1New; y(j,1)=y1New; z(j,1)=z1New x(j,2)=x2New; y(j,2)=y2New; z(j,2)=z2New x(j,3)=x3New; y(j,3)=y3New; z(j,3)=z3New x(j,4)=x4New; y(j,4)=y4New; z(j,4)=z4New x(j,5)=x5New; y(j,5)=y5New; z(j,5)=z5New x(j,6)=x6New; y(j,6)=y6New; z(j,6)=z6New x(j,7)=x7New; y(j,7)=y7New; z(j,7)=z7New x(j,8)=x8New; y(j,8)=y8New; z(j,8)=z8New x(j,9)=x9New; y(j,9)=y9New; z(j,9)=z9New x(j,10)=x10New; y(j,10)=y10New; z(j,10)=z10New endif enddo ! molecule loop !--- Adaptive step size --- if (mod(i,updateInter).eq.0) then ratio = acceptCount /(dble(n)*updateInter) rot_ratio = rot_acceptCount/(dble(n)*updateInter) if (ratio.gt.0.5D0) then maxDisplacement=min(maxDisplacement*1.05D0, box*0.25D0) else maxDisplacement=max(maxDisplacement*0.95D0, 0.01D0*sigma) endif if (rot_ratio.gt.0.5D0) then dPhiMax=min(dPhiMax*1.05D0, pi*0.5D0) else dPhiMax=max(dPhiMax*0.95D0, 0.01D0) endif acceptCount=0.0D0; rot_acceptCount=0.0D0 ! Status line printed only every 1000 steps to keep the log compact ! (step-size tuning above still runs every updateInter steps). if (mod(i,1000).eq.0) then write(*,'(A,I8,A,F8.5,A,F6.3,A,F6.3)') & " Step:",i," MaxDisp:",maxDisplacement, & " AccRat:",ratio," RotRat:",rot_ratio endif endif !--- Energy drift correction --- if (mod(i,1000).eq.0) then call potentialCal(box,sigma,rMin,temp,overlap,vRecalc) if (overlap) write(*,*) "WARNING: Overlap at step",i if (abs(vRecalc-v) > 0.01D0) then write(*,*) "Energy drift at step",i,": correcting; diff=",abs(vRecalc-v) v=vRecalc endif endif !--- Flexibility analysis --- if (mod(i,flexibility_output_freq).eq.0 .and. i>nStep/equilFrac) then call analyze_molecular_flexibility(i, box, trim(outputFlex)) bend_samples=bend_samples+1 endif !--- XYZ output every xyzFreq steps --- if (mod(i,xyzFreq).eq.0) then acceptanceRatio =(acceptance /dble(n))*100.0D0 rot_acceptanceRatio=(rot_acceptance/dble(n))*100.0D0 open(unit=outputUnit, file=trim(outputDat), position='append') open(unit=outputUnit+1, file=trim(outputXyz), position='append') write(outputUnit+1,*) n*atomNo write(outputUnit+1,'(A,7(1X,F12.5))') & "box,temp,step,accR,rotAccR,maxDisp,dPhiMax=", & box,temp,dble(i),acceptanceRatio,rot_acceptanceRatio,maxDisplacement,dPhiMax do w1=1,n do v1=1,atomNo write(outputUnit,*) x(w1,v1),y(w1,v1),z(w1,v1) if (v1.eq.1) then; write(outputUnit+1,*) x(w1,v1),y(w1,v1),z(w1,v1),1,0,0 else if (v1.eq.2) then; write(outputUnit+1,*) x(w1,v1),y(w1,v1),z(w1,v1),0,1,0 else if (v1.eq.9) then; write(outputUnit+1,*) x(w1,v1),y(w1,v1),z(w1,v1),1,0,1 else if (v1.eq.10) then; write(outputUnit+1,*) x(w1,v1),y(w1,v1),z(w1,v1),1,1,0 else ; write(outputUnit+1,*) x(w1,v1),y(w1,v1),z(w1,v1),0,0,1 endif enddo enddo close(outputUnit); close(outputUnit+1) acceptance=0.0D0; rot_acceptance=0.0D0 endif enddo ! MC loop !======================================================================= ! FINAL FLEXIBILITY SUMMARY !======================================================================= if (bend_samples > 0) then avg_end_to_end = sum_end_to_end/dble(bend_samples) avg_rg = sum_rg /dble(bend_samples) write(*,*) "=== FLEXIBILITY SUMMARY ===" write(*,*) "Samples :", bend_samples write(*,*) "Avg Ree :", avg_end_to_end write(*,*) "Avg Rg :", avg_rg write(*,*) "Ree fluctuation :", sqrt(sum_end_to_end_sq/bend_samples - avg_end_to_end**2) write(*,*) "Rg fluctuation :", sqrt(sum_rg_sq/bend_samples - avg_rg**2) write(*,*) "Rg/Ree :", avg_rg/avg_end_to_end write(*,*) "Bond angles (deg):" do i=1,8 avg_bond_angle=sum_bond_angles(i)/dble(bend_samples) write(*,*) " Angle",i,"=",avg_bond_angle*180.0D0/pi enddo persistence_length=bond_length_tail_tail/(1.0D0-cos(sum_bond_angles(1)/dble(bend_samples))) write(*,*) "Persistence length:", persistence_length endif !======================================================================= ! FINAL ENERGY CHECK !======================================================================= call potentialCal(box,sigma,rMin,temp,overlap,vEnd) if (overlap) then; write(*,*) "ERROR: Overlap in final configuration"; stop; endif write(*,*) "Temperature:", temp, " Box:", box if (abs(vEnd-v) .gt. 0.001D0) then write(*,*) "Energy conservation WARNING: running v=",v," vEnd=",vEnd write(*,*) "Relative error =", abs(vEnd-v)/abs(vEnd)*100.0D0,"%" else write(*,*) "Energy conservation: GOOD vEnd=",vEnd endif flush(6) ! force the final summary/energy block out of the stdout buffer end !======================================================================= ! PERIODIC BOUNDARY CONDITIONS !======================================================================= function apply_pbc(coord, box_size) implicit none double precision coord,box_size,apply_pbc apply_pbc = coord - box_size*floor(coord/box_size) return end function apply_pbc !======================================================================= ! FULL SYSTEM POTENTIAL !======================================================================= subroutine potentialCal(box,sigma,rMin,temp,overlap,v) common /position1/x,y,z,ex,ey,ez,n,atomNo common /potparams/ wF,wCut,attraction,frac,w11_param,kappa_param double precision rCut,sigma,rMin,v,x(2000,10),y(2000,10),z(2000,10),xij,yij,zij,rC double precision rijSquare,vIJ,x1,y1,z1,ex(2000,10),ey(2000,10),ez(2000,10),box,boxInverse double precision bTT,wCut,wF,vFlat,temp,frac,attraction,w11_param,kappa_param logical overlap integer n,i,j,k,l,atomNo overlap=.false. boxInverse=1.0D0/box rC=0.30D0*box bTT=1.0D0*sigma rCut=(2.0D0**(1.0D0/6.0D0))*sigma v=0.0D0 do i=1,n do j=i,n do k=1,atomNo x1=x(i,k); y1=y(i,k); z1=z(i,k) do l=1,atomNo if (i/=j .or. (i==j .and. k/=l)) then xij=x1-x(j,l); yij=y1-y(j,l); zij=z1-z(j,l) xij=xij-(anint(xij*boxInverse)*box) yij=yij-(anint(yij*boxInverse)*box) zij=zij-(anint(zij*boxInverse)*box) rijSquare=(xij**2)+(yij**2)+(zij**2) if (rijSquare .le. rMin**2) then; overlap=.true.; return; endif if (sqrt(rijSquare) .le. rC) then vIJ=frac*((w11_param/temp)*exp(-(kappa_param*sqrt(rijSquare))))/(sqrt(rijSquare)) vFlat=0.0D0 if (i/=j) then if (sqrt(rijSquare) .lt. (rCut+wF)) then vFlat=-attraction else if (sqrt(rijSquare).ge.rCut .and. sqrt(rijSquare).le.(wF+wCut)) then vFlat=4.0D0*(((bTT/(sqrt(rijSquare)-wF))**12) & -((bTT/(sqrt(rijSquare)-wF))**6)) endif endif v=v+vIJ+vFlat endif endif enddo enddo enddo enddo return end !======================================================================= ! SINGLE MOLECULE ENERGY !======================================================================= subroutine energyCal(x1Old,y1Old,z1Old,x2Old,y2Old,z2Old,x3Old,y3Old,z3Old, & x4Old,y4Old,z4Old,x5Old,y5Old,z5Old,x6Old,y6Old,z6Old, & x7Old,y7Old,z7Old,x8Old,y8Old,z8Old,x9Old,y9Old,z9Old, & x10Old,y10Old,z10Old,j,box,sigma,temp,v) common /position1/x,y,z,ex,ey,ez,n,atomNo common /potparams/ wF,wCut,attraction,frac,w11_param,kappa_param double precision rCut,sigma,v,xij,yij,zij,rC,rMin double precision rijSquare,x(2000,10),y(2000,10),z(2000,10),x1,y1,z1,x2,y2,z2 double precision xtrial(2000,10),ytrial(2000,10),ztrial(2000,10) double precision ex(2000,10),ey(2000,10),ez(2000,10) double precision bTT,wCut,wF,vFlat,box,boxInverse,vIJ double precision w11_param,kappa_param,frac,attraction,temp double precision x1Old,y1Old,z1Old,x2Old,y2Old,z2Old,x3Old,y3Old,z3Old, & x4Old,y4Old,z4Old,x5Old,y5Old,z5Old,x6Old,y6Old,z6Old, & x7Old,y7Old,z7Old,x8Old,y8Old,z8Old,x9Old,y9Old,z9Old,x10Old,y10Old,z10Old integer j,n,i,k,l,atomNo save xtrial,ytrial,ztrial boxInverse=1.0D0/box; rC=0.3D0*box; rMin=0.7D0*sigma bTT=1.0D0*sigma; rCut=(2.0D0**(1.0D0/6.0D0))*sigma xtrial(j,1)=x1Old; ytrial(j,1)=y1Old; ztrial(j,1)=z1Old xtrial(j,2)=x2Old; ytrial(j,2)=y2Old; ztrial(j,2)=z2Old xtrial(j,3)=x3Old; ytrial(j,3)=y3Old; ztrial(j,3)=z3Old xtrial(j,4)=x4Old; ytrial(j,4)=y4Old; ztrial(j,4)=z4Old xtrial(j,5)=x5Old; ytrial(j,5)=y5Old; ztrial(j,5)=z5Old xtrial(j,6)=x6Old; ytrial(j,6)=y6Old; ztrial(j,6)=z6Old xtrial(j,7)=x7Old; ytrial(j,7)=y7Old; ztrial(j,7)=z7Old xtrial(j,8)=x8Old; ytrial(j,8)=y8Old; ztrial(j,8)=z8Old xtrial(j,9)=x9Old; ytrial(j,9)=y9Old; ztrial(j,9)=z9Old xtrial(j,10)=x10Old; ytrial(j,10)=y10Old; ztrial(j,10)=z10Old v=0.0D0 do i=1,n do k=1,atomNo x1=xtrial(j,k); y1=ytrial(j,k); z1=ztrial(j,k) do l=1,atomNo if (i/=j .or. (i==j .and. k/=l)) then ! For an intramolecular pair use molecule j's TRIAL coordinates ! for both beads. Otherwise a proposed (not-yet-committed) move ! mixes trial bead k with the old global bead l, producing a ! spurious intramolecular energy difference that accumulates as ! a large drift in the running energy v. if (i==j) then x2=xtrial(j,l); y2=ytrial(j,l); z2=ztrial(j,l) else x2=x(i,l); y2=y(i,l); z2=z(i,l) endif xij=x1-x2; yij=y1-y2; zij=z1-z2 xij=xij-(anint(xij*boxInverse)*box) yij=yij-(anint(yij*boxInverse)*box) zij=zij-(anint(zij*boxInverse)*box) rijSquare=(xij**2)+(yij**2)+(zij**2) if (rijSquare .le. rMin**2) then; v=1.0D15; return; endif if (sqrt(rijSquare) .le. rC) then vIJ=frac*((w11_param/temp)*exp(-(kappa_param*sqrt(rijSquare))))/(sqrt(rijSquare)) vFlat=0.0D0 if (i/=j) then if (sqrt(rijSquare) .lt. (rCut+wF)) then vFlat=-attraction else if (sqrt(rijSquare).ge.rCut .and. sqrt(rijSquare).le.(wF+wCut)) then vFlat=4.0D0*(((bTT/(sqrt(rijSquare)-wF))**12) & -((bTT/(sqrt(rijSquare)-wF))**6)) endif endif v=v+vIJ+vFlat endif endif enddo enddo enddo return end !======================================================================= ! OVERLAP TEST !======================================================================= subroutine test(x1New,y1New,z1New,x2New,y2New,z2New,x3New,y3New,z3New, & x4New,y4New,z4New,x5New,y5New,z5New,x6New,y6New,z6New, & x7New,y7New,z7New,x8New,y8New,z8New,x9New,y9New,z9New, & x10New,y10New,z10New,j,box,rMin,overlap) common /position1/x,y,z,ex,ey,ez,n,atomNo double precision rMin,x(2000,10),y(2000,10),z(2000,10),xij,yij,zij,x1,y1,z1 double precision rijSquare,ex(2000,10),ey(2000,10),ez(2000,10) double precision xtrial(2000,10),ytrial(2000,10),ztrial(2000,10),box,boxInverse double precision x1New,y1New,z1New,x2New,y2New,z2New,x3New,y3New,z3New, & x4New,y4New,z4New,x5New,y5New,z5New,x6New,y6New,z6New, & x7New,y7New,z7New,x8New,y8New,z8New,x9New,y9New,z9New,x10New,y10New,z10New logical overlap integer n,i,j,k,l,atomNo save xtrial,ytrial,ztrial overlap=.false.; boxInverse=1.0D0/box xtrial(j,1)=x1New; ytrial(j,1)=y1New; ztrial(j,1)=z1New xtrial(j,2)=x2New; ytrial(j,2)=y2New; ztrial(j,2)=z2New xtrial(j,3)=x3New; ytrial(j,3)=y3New; ztrial(j,3)=z3New xtrial(j,4)=x4New; ytrial(j,4)=y4New; ztrial(j,4)=z4New xtrial(j,5)=x5New; ytrial(j,5)=y5New; ztrial(j,5)=z5New xtrial(j,6)=x6New; ytrial(j,6)=y6New; ztrial(j,6)=z6New xtrial(j,7)=x7New; ytrial(j,7)=y7New; ztrial(j,7)=z7New xtrial(j,8)=x8New; ytrial(j,8)=y8New; ztrial(j,8)=z8New xtrial(j,9)=x9New; ytrial(j,9)=y9New; ztrial(j,9)=z9New xtrial(j,10)=x10New; ytrial(j,10)=y10New; ztrial(j,10)=z10New do i=1,n if (i/=j) then do k=1,atomNo x1=xtrial(j,k); y1=ytrial(j,k); z1=ztrial(j,k) do l=1,atomNo xij=x1-x(i,l); yij=y1-y(i,l); zij=z1-z(i,l) xij=xij-(anint(xij*boxInverse)*box) yij=yij-(anint(yij*boxInverse)*box) zij=zij-(anint(zij*boxInverse)*box) if ((xij**2+yij**2+zij**2).le.rMin**2) then; overlap=.true.; return; endif enddo enddo endif enddo return end !======================================================================= ! FLEXIBLE ROTATION (pivot move) !======================================================================= subroutine rotate_flexible_molecule(j, box, rMin, bond_length, max_attempts) implicit none integer j,max_attempts,n,atomNo double precision box,rMin,bond_length common /position1/x,y,z,ex,ey,ez,n,atomNo double precision ex(2000,10),ey(2000,10),ez(2000,10),x(2000,10),y(2000,10),z(2000,10) double precision pi,boxInverse,apply_pbc,theta,phi,new_x,new_y,new_z,old_x,old_y,old_z double precision bvx,bvy,bvz integer k,attempt logical overlap,bond_ok pi=4.0D0*atan(1.0D0); boxInverse=1.0D0/box do k=2,10 old_x=x(j,k); old_y=y(j,k); old_z=z(j,k) do attempt=1,max_attempts call random_number(theta); call random_number(phi) theta=theta*2.0D0*pi; phi=phi*pi new_x=apply_pbc(x(j,k-1)+bond_length*sin(phi)*cos(theta),box) new_y=apply_pbc(y(j,k-1)+bond_length*sin(phi)*sin(theta),box) new_z=apply_pbc(z(j,k-1)+bond_length*cos(phi),box) call check_single_bead_overlap(j,k,new_x,new_y,new_z,box,rMin,overlap) bvx=new_x-x(j,k-1); bvy=new_y-y(j,k-1); bvz=new_z-z(j,k-1) bvx=bvx-box*anint(bvx*boxInverse) bvy=bvy-box*anint(bvy*boxInverse) bvz=bvz-box*anint(bvz*boxInverse) bond_ok=(abs(sqrt(bvx**2+bvy**2+bvz**2)-bond_length)<0.01D0) if (.not.overlap .and. bond_ok) then x(j,k)=new_x; y(j,k)=new_y; z(j,k)=new_z; exit endif enddo if (attempt>max_attempts) then x(j,k)=old_x; y(j,k)=old_y; z(j,k)=old_z endif enddo end subroutine rotate_flexible_molecule !======================================================================= ! SINGLE BEAD OVERLAP CHECK !======================================================================= subroutine check_single_bead_overlap(j,k,test_x,test_y,test_z,box,rMin,overlap) implicit none integer j,k,n,atomNo double precision test_x,test_y,test_z,box,rMin logical overlap common /position1/x,y,z,ex,ey,ez,n,atomNo double precision ex(2000,10),ey(2000,10),ez(2000,10),x(2000,10),y(2000,10),z(2000,10) integer i,l double precision xij,yij,zij,rij_sq,boxInverse overlap=.false.; boxInverse=1.0D0/box do i=1,n if (i/=j) then do l=1,atomNo xij=test_x-x(i,l); yij=test_y-y(i,l); zij=test_z-z(i,l) xij=xij-box*anint(xij*boxInverse) yij=yij-box*anint(yij*boxInverse) zij=zij-box*anint(zij*boxInverse) if ((xij**2+yij**2+zij**2)<=rMin**2) then; overlap=.true.; return; endif enddo endif enddo do l=1,atomNo if (l/=k .and. abs(l-k)>1) then xij=test_x-x(j,l); yij=test_y-y(j,l); zij=test_z-z(j,l) xij=xij-box*anint(xij*boxInverse) yij=yij-box*anint(yij*boxInverse) zij=zij-box*anint(zij*boxInverse) if ((xij**2+yij**2+zij**2)<=rMin**2) then; overlap=.true.; return; endif endif enddo end subroutine check_single_bead_overlap !======================================================================= ! FLEXIBILITY ANALYSIS !======================================================================= subroutine analyze_molecular_flexibility(step, box, flexFile) implicit none integer step,n,atomNo double precision box character*(*) flexFile common /position1/x,y,z,ex,ey,ez,n,atomNo double precision ex(2000,10),ey(2000,10),ez(2000,10),x(2000,10),y(2000,10),z(2000,10) double precision end_to_end_dist,radius_gyration,bond_angles(8) double precision sum_end_to_end,sum_rg,sum_end_to_end_sq,sum_rg_sq double precision avg_bond_angle,sum_bond_angles(9) integer bend_samples common /flexibility_stats/ sum_end_to_end,sum_rg,sum_end_to_end_sq,sum_rg_sq,sum_bond_angles,bend_samples integer i,j double precision dx,dy,dz,dist,boxInverse,cx,cy,cz double precision rg_x,rg_y,rg_z,v1x,v1y,v1z,v2x,v2y,v2z double precision dot_prod,angle,mag1,mag2,pi,curv,tot_curv,cont_len,avg_ree,avg_rg_mol pi=4.0D0*atan(1.0D0); boxInverse=1.0D0/box avg_ree=0.0D0; avg_rg_mol=0.0D0 open(unit=200, file=flexFile, position='append') do i=1,n dx=x(i,10)-x(i,1); dy=y(i,10)-y(i,1); dz=z(i,10)-z(i,1) dx=dx-box*anint(dx*boxInverse); dy=dy-box*anint(dy*boxInverse); dz=dz-box*anint(dz*boxInverse) end_to_end_dist=sqrt(dx**2+dy**2+dz**2) avg_ree=avg_ree+end_to_end_dist cx=0.0D0; cy=0.0D0; cz=0.0D0 do j=1,10; cx=cx+x(i,j); cy=cy+y(i,j); cz=cz+z(i,j); enddo cx=cx/10.0D0; cy=cy/10.0D0; cz=cz/10.0D0 rg_x=0.0D0; rg_y=0.0D0; rg_z=0.0D0 do j=1,10 dx=x(i,j)-cx; dy=y(i,j)-cy; dz=z(i,j)-cz dx=dx-box*anint(dx*boxInverse); dy=dy-box*anint(dy*boxInverse); dz=dz-box*anint(dz*boxInverse) rg_x=rg_x+dx**2; rg_y=rg_y+dy**2; rg_z=rg_z+dz**2 enddo radius_gyration=sqrt((rg_x+rg_y+rg_z)/10.0D0) avg_rg_mol=avg_rg_mol+radius_gyration do j=2,9 v1x=x(i,j)-x(i,j-1); v1y=y(i,j)-y(i,j-1); v1z=z(i,j)-z(i,j-1) v2x=x(i,j+1)-x(i,j); v2y=y(i,j+1)-y(i,j); v2z=z(i,j+1)-z(i,j) v1x=v1x-box*anint(v1x*boxInverse); v1y=v1y-box*anint(v1y*boxInverse); v1z=v1z-box*anint(v1z*boxInverse) v2x=v2x-box*anint(v2x*boxInverse); v2y=v2y-box*anint(v2y*boxInverse); v2z=v2z-box*anint(v2z*boxInverse) dot_prod=v1x*v2x+v1y*v2y+v1z*v2z mag1=sqrt(v1x**2+v1y**2+v1z**2); mag2=sqrt(v2x**2+v2y**2+v2z**2) if (mag1>0.0D0 .and. mag2>0.0D0) then bond_angles(j-1)=acos(max(-1.0D0,min(1.0D0,dot_prod/(mag1*mag2)))) else bond_angles(j-1)=0.0D0 endif enddo tot_curv=0.0D0; cont_len=0.0D0 do j=1,9 dx=x(i,j+1)-x(i,j); dy=y(i,j+1)-y(i,j); dz=z(i,j+1)-z(i,j) dx=dx-box*anint(dx*boxInverse); dy=dy-box*anint(dy*boxInverse); dz=dz-box*anint(dz*boxInverse) dist=sqrt(dx**2+dy**2+dz**2); cont_len=cont_len+dist if (j>1) then; curv=abs(pi-bond_angles(j-1))/dist; tot_curv=tot_curv+curv; endif enddo write(200,'(I8,I5,6F12.6)') step,i,end_to_end_dist,radius_gyration,cont_len,tot_curv, & end_to_end_dist/cont_len, radius_gyration/end_to_end_dist enddo avg_ree=avg_ree/n; avg_rg_mol=avg_rg_mol/n sum_end_to_end=sum_end_to_end+avg_ree; sum_rg=sum_rg+avg_rg_mol sum_end_to_end_sq=sum_end_to_end_sq+avg_ree**2; sum_rg_sq=sum_rg_sq+avg_rg_mol**2 do j=1,8 avg_bond_angle=0.0D0 do i=1,n v1x=x(i,j+1)-x(i,j); v1y=y(i,j+1)-y(i,j); v1z=z(i,j+1)-z(i,j) v2x=x(i,j+2)-x(i,j+1); v2y=y(i,j+2)-y(i,j+1); v2z=z(i,j+2)-z(i,j+1) v1x=v1x-box*anint(v1x*boxInverse); v1y=v1y-box*anint(v1y*boxInverse); v1z=v1z-box*anint(v1z*boxInverse) v2x=v2x-box*anint(v2x*boxInverse); v2y=v2y-box*anint(v2y*boxInverse); v2z=v2z-box*anint(v2z*boxInverse) dot_prod=v1x*v2x+v1y*v2y+v1z*v2z mag1=sqrt(v1x**2+v1y**2+v1z**2); mag2=sqrt(v2x**2+v2y**2+v2z**2) if (mag1>0.0D0 .and. mag2>0.0D0) then angle=acos(max(-1.0D0,min(1.0D0,dot_prod/(mag1*mag2)))) avg_bond_angle=avg_bond_angle+angle endif enddo sum_bond_angles(j)=sum_bond_angles(j)+(avg_bond_angle/n) enddo write(*,'(A,I8,A,F8.4,A,F8.4)') "Step",step,": Avg Ree=",avg_ree,", Avg Rg=",avg_rg_mol close(200) end subroutine analyze_molecular_flexibility