MODULE global_variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Contains all the global variables that will be used in the program ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Number of particles !----------------------------------------------------------------------! ! Number of nitrogen molecules INTEGER :: n_molecules INTEGER,PARAMETER :: n_molecules_max=8000 ! Positions !----------------------------------------------------------------------! ! Positions of the centers of mass of the N molec. REAL,DIMENSION(n_molecules_max) :: rx, ry, rz ! Geometry !----------------------------------------------------------------------! ! Simulation box length REAL :: box ! Physical variables !----------------------------------------------------------------------! ! Temperature (non-dimensional) REAL :: temp ! Energy, etc REAL :: energy ! Block accumulators REAL :: blkener, blkpress ! Virial and pressure REAL :: pressure ! Parameters !----------------------------------------------------------------------! ! Cutoff distances REAL :: cutoff_nn, cutoff_nn_square ! Very small distances will cause huge energies REAL, PARAMETER :: huge_number = 10000 REAL :: small_rsq_nn, huge_energy ! Acceptance/Rejections Statistics !----------------------------------------------------------------------! ! Translation statistics INTEGER :: acctrans, ntrans ! Monte Carlo parameters !----------------------------------------------------------------------! ! Number of blocks INTEGER :: nblocks ! Number of steps per block INTEGER :: nsteps ! Maximum displacement in a MC move REAL :: maxdisp ! Adjust the parameter maxdisp every nadjust blocks INTEGER :: nadjust ! Radial distribution function !----------------------------------------------------------------------! INTEGER,PARAMETER :: maxbin = 100 INTEGER, DIMENSION(maxbin) :: hist REAL :: delr,constant INTEGER :: nsamplegr,totalsamplesgr ! Others !----------------------------------------------------------------------! ! Dummy variable for random number generator INTEGER :: idum ! Random number generator REAL,EXTERNAL :: random ! Pi REAL, PARAMETER :: pi=3.14159 END MODULE global_variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PROGRAM nvt !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Main program ! ! Procedures used: ! ! - start, runner, finish ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE CALL start CALL runner CALL finish STOP END PROGRAM nvt !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE runner !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Performs nblocks blocks of nsteps Monte Carlo steps ! ! ! ! Procedures used: ! ! - random ! ! - translate ! ! - total_energy ! ! - adjust_mc_parameters ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Local variables INTEGER :: iblock, istep ! Calculate the total energy CALL total_energy(energy,pressure) ! Perform runs over blocks DO iblock=1,nblocks ! Set block accumulators to zero blkener = 0. blkpress = 0. ! Perform runs over steps DO istep=1,nsteps ! Perform an MC step CALL translate ! Update block accumulators blkener = blkener + energy blkpress = blkpress + pressure ! End loop over steps END DO ! Average properties at the end of the block blkener = blkener / REAL(nsteps) blkpress = blkpress / REAL(nsteps) ! Write results in output file WRITE(40,FMT='(I5,5X,3F18.6)') & &iblock,blkener/n_molecules,energy/n_molecules,& &(n_molecules*temp-blkpress)/box**3 ! Adjust MC parameters IF (MOD(iblock,nadjust) == 0) CALL adjust_mc_parameters ! Sample radial distribution function IF (MOD(iblock,nsamplegr) == 0) CALL sample_gr ! End loop over blocks END DO ! Close output file CLOSE(UNIT=40) RETURN END SUBROUTINE runner !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! REAL FUNCTION random(idum) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Generates a random number between 0.0 and 1.0 ! ! Modification of ran3, taken from "Numerical Recipes in Fortran" ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE ! Local variables REAL, PARAMETER :: mbig=4000000., mseed=1618033., mz=0., fac=1./mbig INTEGER :: idum,i, ii, k INTEGER, SAVE :: iff, inext, inextp REAL :: mj, mk REAL, SAVE :: ma(55) DATA iff /0/ ! Initialization IF(idum<0.OR.iff==0) THEN iff=1 mj=mseed-IABS(idum) mj=MOD(mj, mbig) ma(55)=mj mk=1 DO i=1,54 ii=MOD(21*i,55) ma(ii)=mk mk=mj-mk IF(mk= (huge_energy-1.)) THEN accepted = .FALSE. ELSEIF (random(idum) < EXP(-deltav/temp)) THEN accepted = .TRUE. ELSE accepted = .FALSE. END IF ! If the move is accepted... IF (accepted) THEN ! Update the number of accepted translation moves acctrans = acctrans + 1 ! Update positions rx(imol) = rx_new ry(imol) = ry_new rz(imol) = rz_new ! Update the value of the energy energy = energy + deltav ! Update the value of pressure pressure = pressure + virialnew - virialold END IF ! Update the number of translation moves ntrans = ntrans + 1 RETURN END SUBROUTINE translate !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE energy_calculation(itest,rxi,ryi,rzi,potenergy,virial) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Calculates the potential energy of molecule"i" (test molecule) ! ! with all the other molecules/atoms ("j") ! ! ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Variables with argument association INTEGER :: itest REAL :: rxi,ryi,rzi,potenergy,virial ! Local variables REAL :: energy_ij,virial_ij REAL :: factor REAL :: rxij,ryij,rzij,rijsq INTEGER :: j ! Contribution due to interactions with other nitrogen molecules !Initialize variables potenergy = 0. virial = 0. ! Start loop over all nitrogen molecules DO j=1,n_molecules ! avoid self-interactions IF (j==itest) CYCLE ! calculate vector between centers of mass rxij = rx(j) - rxi ryij = ry(j) - ryi rzij = rz(j) - rzi ! apply minimum image convention rxij = rxij - ANINT(rxij/box)*box ryij = ryij - ANINT(ryij/box)*box rzij = rzij - ANINT(rzij/box)*box ! Calculate the distance between centers of mass rijsq = rxij**2 + ryij**2 + rzij**2 ! Apply cutoff distance IF (rijsq <= cutoff_nn_square) THEN ! Calculate the potential due to dispersion interactions ! Check if the distance is very small IF (rijsq <= small_rsq_nn) THEN potenergy = huge_energy RETURN END IF ! Evaluate the L-J potential factor = (1./ rijsq)**3 energy_ij = factor*(factor-1) ! Evaluate the virial virial_ij = factor*(1-2*factor) ! Add contributions to energy and virial potenergy = potenergy + energy_ij virial = virial + virial_ij ! End the calculations for distance smaller than cutoff END IF ! End loop over nitrogen molecules END DO ! Multiply energy and virial by the L-J factor potenergy = 4. * potenergy virial = 8. * virial END SUBROUTINE energy_calculation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE total_energy(potenergy,virial) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Calculates the total energy of the system ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Variables with argument association REAL :: potenergy,virial ! Local variables REAL :: rxi,ryi,rzi REAL :: energy_ij, virial_ij REAL :: factor REAL :: rxij,ryij,rzij,rijsq INTEGER :: j,imol ! Initialize parameters potenergy = 0. virial = 0. ! Start loop over all nitrogen molecules DO imol=1,n_molecules-1 rxi = rx(imol) ryi = ry(imol) rzi = rz(imol) ! Contribution due to interactions with other nitrogen molecules ! Start loop over other nitrogen molecules DO j=imol+1,n_molecules ! calculate vector between centers of mass rxij = rx(j) - rxi ryij = ry(j) - ryi rzij = rz(j) - rzi ! apply minimum image convention rxij = rxij - ANINT(rxij/box)*box ryij = ryij - ANINT(ryij/box)*box rzij = rzij - ANINT(rzij/box)*box ! Calculate the distance between centers of mass rijsq = rxij**2 + ryij**2 + rzij**2 ! Apply cutoff distance IF (rijsq <= cutoff_nn_square) THEN ! Calculate the potential due to dispersion interactions ! Check if the distance is very small IF (rijsq <= small_rsq_nn) THEN energy = huge_energy PRINT*, 'ERROR: Try with a different initial config' STOP END IF ! Evaluate the L-J potential factor = (1./ rijsq)**3 energy_ij = factor*(factor-1) ! Calculate the virial virial_ij = factor*(1-2*factor) ! Add contributions to energy and pressure potenergy = potenergy + energy_ij virial = virial + virial_ij ! End the calculations for distance smaller than cutoff END IF ! End loop over nitrogen molecules END DO ! End loop over all nitrogen molecules END DO ! Multiply energy and pressure by the L-J factor potenergy = 4. * potenergy virial = 8. * virial END SUBROUTINE total_energy !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE adjust_mc_parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Adjusts MC parameter maxdisp ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Adjust maxdisp IF ( (REAL(acctrans)/REAL(ntrans)) > 0.5) THEN IF (maxdisp < box) maxdisp = 1.05 * maxdisp ELSE maxdisp = 0.95 * maxdisp ENDIF RETURN END SUBROUTINE adjust_mc_parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE sample_gr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Updates the histogram for g(r) ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Local variables INTEGER :: imol,j REAL :: rxi,ryi,rzi REAL :: rxij,ryij,rzij,rijsq,rij INTEGER :: bin totalsamplesgr = totalsamplesgr + 1 ! Start loop over all nitrogen molecules DO imol=1,n_molecules-1 rxi = rx(imol) ryi = ry(imol) rzi = rz(imol) ! Contribution due to interactions with other nitrogen molecules ! Start loop over other molecules DO j=imol+1,n_molecules ! calculate vector between centers of mass rxij = rx(j) - rxi ryij = ry(j) - ryi rzij = rz(j) - rzi ! apply minimum image convention rxij = rxij - ANINT(rxij/box)*box ryij = ryij - ANINT(ryij/box)*box rzij = rzij - ANINT(rzij/box)*box ! Calculate the distance between centers of mass rijsq = rxij**2 + ryij**2 + rzij**2 rij = SQRT(rijsq) IF (rij < box/2) THEN bin = INT(rij/delr) + 1 hist(bin) = hist(bin) + 2 END IF ! End loop over other molecules END DO ! End loop over molecules END DO RETURN END SUBROUTINE sample_gr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE start !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Sets everything up to start running the MC code ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Local variables INTEGER :: imol,bin ! Read the input file OPEN(UNIT=10,FILE='inputfile',STATUS='OLD',ACCESS='SEQUENTIAL',& &ACTION='READ') READ(10,*) READ(10,*) temp READ(10,*) cutoff_nn cutoff_nn_square = cutoff_nn * cutoff_nn READ(10,*) READ(10,*) READ(10,*) nblocks READ(10,*) nsteps READ(10,*) nadjust READ(10,*) nsamplegr READ(10,*) maxdisp READ(10,*) box CLOSE(UNIT=10) ! Calculate the very small distances that will cause huge energies huge_energy = temp * huge_number small_rsq_nn = (0.8)**2 ! Get initial positions of nitrogen molecules OPEN(UNIT=20,FILE='initial.conf',STATUS='OLD',ACCESS='SEQUENTIAL',& &ACTION='READ') READ(20,*) n_molecules READ(20,*) READ(20,*) READ(20,*) ! Read positions DO imol=1,n_molecules READ(20,*) rx(imol),ry(imol),rz(imol) END DO CLOSE(UNIT=20) ! Initialize dummy variable for random number generator CALL initialize_idum(idum) ! Initialize acc/rej statistics acctrans = 0 ntrans = 0 ! Initialize radial distribution function DO bin=1,maxbin hist(bin) = 0 END DO delr = box/(2*maxbin) constant=4*pi*n_molecules/(3*box**3) totalsamplesgr = 0 ! Open file to write results OPEN(UNIT=40,FILE='blkaverages',STATUS='UNKNOWN',ACCESS='SEQUENTIAL') WRITE(40,*) 'NVT Results (blocks averages)' WRITE(40,*) WRITE(40,FMT='(F14.6, '' = Temperature'')') temp WRITE(40,FMT='(9X,I5, '' = Number of molecules'')') n_molecules WRITE(40,FMT='(F14.6, '' = Volume of the simulation box'')') box**3 WRITE(40,*) WRITE(40,*) 'block Ublock Uinst P' RETURN END SUBROUTINE start !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE finish !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Reports the final results of the MC code ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Jorge Pikunic ! ! North Carolina State University ! ! Department of Chemical Engineering ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE global_variables IMPLICIT NONE ! Local variables INTEGER :: imol,j REAL :: rlow,rup,nideal INTEGER :: bin OPEN(UNIT=50,FILE='final.conf',STATUS='UNKNOWN',ACCESS='SEQUENTIAL') WRITE(50,FMT='(1I4,12X,''Number of nitrogen molecules'')') n_molecules WRITE(50,*) ' ' WRITE(50,*) ' POSITION OF EACH MOLECULE' WRITE(50,*) ' rx ry rz ' DO imol=1,n_molecules WRITE(50,FMT='(3F12.6)') rx(imol),ry(imol),rz(imol) END DO CLOSE(UNIT=50) ! Calculate radial distribution function OPEN(UNIT=60,FILE='gr',STATUS='UNKNOWN',ACCESS='SEQUENTIAL') DO bin=1,maxbin rlow = delr*(bin-1) rup = rlow + delr nideal = constant*(rup**3-rlow**3) WRITE(60,FMT='(F5.3,F12.6)') rlow + delr/2 , & &hist(bin)/nideal/REAL(n_molecules)/REAL(totalsamplesgr) END DO CLOSE(UNIT=60) RETURN END SUBROUTINE finish