PROGRAM AMMONIA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! USER AGREEMENT ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Except as otherwise expressly permitted users of this code may not:! ! (i) modify or create any derivative works of this code or ! ! documentation, including translation or localization; ! ! (ii) redistribute, encumber, sell, rent, lease, sublicense, or ! ! otherwise transfer rights to this code; (iii) remove or alter any ! ! trademark, logo, copyright or other proprietary notices, legends, ! ! symbols or labels in this code; or (iv) publish any results ! ! obtained with this code without the author's prior written consent.! ! Download or usage of this code constitutes agreement with the above! ! restrictions. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! VARIABLES ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! N: total # of phases present in the system ! ! nph: specific phase, bulk (nph=1), slit (nph=2) ! ! xa: atomic positions ! ! xm: molecular center of mass ! ! ex: orientation of LJ cite in a molecule ! ! qx: orientation of partial charge in a molecule ! ! bl: length of COM to LJ site ! ! ql: length of COM to partial charge site ! ! Dgamax_b: maximum attempted rotation (radians) in bulk ! ! Dgamax_s: maximum attempted rotation (radians) in slit ! ! disp_b: maximum attempted displacement (sigma) in bulk ! ! disp_s: maximum attempted displacement (sigma) in slit ! ! t_equil: number of steps before averages are taken ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' CALL Zero CALL Make_file temp_cycle: DO kk = TP_i,TP_f ! cycle through T/P conditions CALL Initialize ! Read the T/P parameters, and initialize IF(n_cont.EQ.1) CALL Continue_Run ! continue from previous configuration IF(kk.EQ.TP_i) CALL Animate ! write xyz file if first run nph=1; CALL Partition nph=2; CALL Partition main_loop: DO WHILE (t.LT.tmax) ! Begin the main loop nTot_mol_1=0 nTot_mol_2=0 count: DO i=1,nspecies ! Count # of mol in each phase nTot_mol_1 = nTot_mol_1 + Nb(i,1) nTot_mol_2 = nTot_mol_2 + Nb(i,2) END DO count zz = REAL(nTot_mol_1)/REAL(nTot_mol_1+nTot_mol_2) ! Choose a phase to work in IF(ran3(iseed).LT.zz) THEN; nph=1; ELSE; nph=2; END IF zzz = ran3(iseed) ! randomly select step IF (zzz.LT.0.30) THEN; CALL Move ! moves a molecule ELSE IF(zzz.LT.0.40) THEN; CALL Rotate ! rotates a molecule ELSE IF(zzz.LT.0.70) THEN ! forward reaction step CALL Forward ELSE IF(zzz.LT.1.00) THEN ! reverse reaction step CALL Reverse END IF CALL Averages IF (loop.EQ.500) loop = 0 IF (loop.EQ.0) THEN CALL Show_status CALL Boundary CALL Volume_NPT ! adjust the volume for an NPT simulation CALL Averages CALL RDF nph=1; CALL Partition nph=2; CALL Partition END IF CALL Chemical_potential ! calculate chemical potential ! zzz = ran3(iseed) ! Exchange particles ! IF(zzz.LT.0.5) CALL Exchange_1_2 ! Move molecule from bulk to pore ! IF(zzz.GE.0.5) CALL Exchange_2_1 ! Move molecule from pore to bulk ! CALL Averages t = t + 1.0 ; loop = loop + 1 ! update cycle # END DO main_loop nph=1 CALL Animate ! write a '.xyz' file nph=1; CALL Partition nph=2; CALL Partition CALL Print_xa ! saves system info. CALL Print_config ! prints configuration file IF (Nrdf.EQ.1) CALL RDF_Print ! Print current RDF data END DO temp_cycle END PROGRAM AMMONIA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Zero INCLUDE 'ammonia.inc' CHARACTER*20 file_name box = 0.0 Na=0; Nb=0; Nc=0; bl=0.0 num_mol=0 acc1=0.0; acc2=0.0; acc3_b=0.0; acc3_s=0.0; acc4_b=0.0; acc4_s=0.0 att1=0.0; att2=0.0; att3_b=0.0; att3_s=0.0; att4_b=0.0; att4_s=0.0 sd1=0.0;sd2=0.0;sd3=0.0;sd4=0.0;sd5=0.0;sd6=0.0;sd7=0.0;sd8=0.0 xa=0.0;ya=0.0;za=0.0;xm=0.0;ym=0.0;zm=0.0;ex=0.0;ey=0.0;ez=0.0 rho_m=0.0; rho_a=0.0;bl=0.0 Dgamax_b=pi/4.0; Dgamax_s=pi/4.0 npress=0; novr=0; P_temp=0.0 sig=0.0;eps=0.0;sig_m=0.0;eps_m=0.0 nins=0; ndel=0; nbias=0 disp_b=1.0; disp_s=1.0; dbox=2.0 cumNb_b = 0.0; cumNb_s = 0.0; numave = 0 cum_B_vir = 0.0 Vol=0.0; P_set=0.0; Temp=0.0; Uold=0.0; Vir=0.0 histC1=0.0; histC2=0.0; histA1=0.0; histA2=0.0; histA3=0.0 histA=0.0;histB=0.0 t = 1.0 loop = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! inert = 0 ! 0 = pure N2,H2,NH3 ; 1 = addition of inert components n_cont = 0 ! 0 = new run, 1 = start from previous configuration ! Check the continue_run subroutine if n_cont=1 box(3,2) = 5.0 ! Slit width: C-C center num_run = 1 ! Output file name Niso = 0 ! Generate an isotherm Nprint = 0 ! Print to screen Nrdf = 0 ! Radial distribution function Nanim = 1 ! Animation file TP_i = 1 ! Temperature # to start with Ntemps = 1 ! Number of temperatures to run Ntemps2 = 1 ! Number of lines read to from file t_equil = 1.0e4 ! # of equilibration moves tmax = 4.0e4 ! total time cut = 4.0 ! Cutoff distance for LJ interactions el_cut = 4.0 ! Cutoff distance for electrostatic interactions V_pore = 2.00 ! Total pore volume num_mol(1)= 2 ! # of molecules of type (1) per side (change this for different size) num_mol(2)= 6 ! # of molecules of type (2) per side num_mol(3)= 0 ! # of molecules of type (3) per side !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! OPEN(unit=120,file="tp_input.dat") ! source file for T/P data DO i=1,Ntemps2 ! read in T/P data READ (120,*) Temp(i) READ (120,*) P_set(i) READ (120,*) Rho_set(i) READ (120,*) t_equil READ (120,*) tmax END DO REWIND(120) PRINT * PRINT 120,'Equilibration time:',t_equil PRINT 120,'Total sim. time: ',tmax num_mol(4)= num_mol(1)+num_mol(2)+num_mol(3) ! length of a side num_mol(5)= num_mol(4)*num_mol(4)*num_mol(4) ! total # of initial molecules nblock = INT((tmax-t_equil)/10.0) ! block length for averages ! PRINT *,'Block length for averages:',nblock ! print block length ! PRINT *,'Run #: ',num_run ! PRINT *,'Pore Width: ',box(3,2) sqcut = cut*cut el_sqcut = el_cut*el_cut TP_f = TP_i + Ntemps - 1 iseed = -2 - num_run ! random number seed D = 251.9 ! change in bond energy 120 FORMAT (1x,A,F12.2) END SUBROUTINE Zero !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Make_file INCLUDE 'ammonia.inc' OPEN(unit=222,file='fnames') ! Generate filenames for run WRITE (222,101) 'am_output',num_run,'.dat' ! unit = 110 WRITE (222,102) 'am_press',num_run,'.dat' ! unit = WRITE (222,103) 'am_status',num_run,'.dat' ! unit = 300 WRITE (222,104) 'am_xa',num_run,'.dat' ! unit = 400 WRITE (222,105) 'am_vol',num_run,'.dat' ! unit = WRITE (222,106) 'xa_am',num_run,'.cycle' ! unit = 600 WRITE (222,107) 'am_grA',num_run,'.dat' ! unit = 700 WRITE (222,107) 'am_grB',num_run,'.dat' ! unit = 800 WRITE (222,103) 'anim_B_am',num_run,'.xyz' ! unit = 900 WRITE (222,103) 'anim_S_am',num_run,'.xyz' ! unit = 1000 WRITE (222,107) 'am_s_E',num_run,'.dat' ! unit = WRITE (222,102) 'am_s_num',num_run,'.dat' ! unit = 1400 WRITE (222,107) 'am_grC',num_run,'.dat' ! unit = 1500 WRITE (222,103) 'am_config',num_run,'.dat' ! unit = 1600 REWIND(222) READ(222,200) fname1; READ(222,200) fname2; READ(222,200) fname3 READ(222,200) fname4; READ(222,200) fname5; READ(222,200) fname6 READ(222,200) fname7; READ(222,200) fname8; READ(222,200) fname9 READ(222,200) fname10; READ(222,200) fname12 READ(222,200) fname14; READ(222,200) fname15;READ(222,200) fname16 OPEN(unit=110,file=fname1) OPEN(unit=300,file=fname3) OPEN(unit=400,file=fname4) ! OPEN(unit=600,file=fname6) ! IF (Nrdf.EQ.1) OPEN(unit=700,file=fname7) ! IF (Nrdf.EQ.1) OPEN(unit=800,file=fname8) ! IF (Nanim.EQ.1) OPEN(unit=900,file=fname9) ! IF (Nanim.EQ.1) OPEN(unit=1000,file=fname10) ! OPEN(unit=1400,file=fname14) ! IF (Nrdf.EQ.1) OPEN(unit=1500,file=fname15) ! OPEN(unit=1600,file=fname16) CLOSE(222); CLOSE(120); REWIND(1100) !!!!!!!!!!!!!!!!!!!!!!!!!! 101 FORMAT(1x,A,I3.3,A4) 102 FORMAT(1x,A,I4.4,A4) 103 FORMAT(1x,A,I3.3,A4) 104 FORMAT(1x,A,I7.7,A4) 105 FORMAT(1x,A,I6.6,A4) 106 FORMAT(1x,A,I5.5,A6) 107 FORMAT(1x,A,I6.6,A4) 108 FORMAT(1x,A,I5.5,A4) 200 FORMAT(1x,A16) END SUBROUTINE Make_file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Initialize INCLUDE 'ammonia.inc' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Model for ammonia is a single LJ sphere with four point charges. A negative ! ! charge is placed on nitrogen and three positive charges are added to ! ! account for the hydrogens. The model is taken from The Journal of Physical ! ! Chemistry, Vol.97, pp. 9241-9247, 1993, by Jiali Gao. ! ! ! ! All of the specified parameters are reduced with the LJ parameters for a ! ! single nitrogen atom : sigma = 3.32 ang and epsilon = 36.4 K ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! The following sections contain the molecular descriptions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! NITROGEN ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! N2_model=1 ! specify the model for N2 to use IF(N2_model.EQ.1) THEN ! USE THIS charge(1,1) = 0.81010 ! partial charge on the bond charge(2,1) = -0.40505 ! partial charge on N #1 charge(3,1) = -0.40505 ! partial charge on N #2 bl(1,1) = 0.16566 ! position of N #1 bl(2,1) = 0.16566 ! position of N #2 ql(1,1) = 0.000 ! charge #1 is located at COM of N2 ql(2,1) = 0.16536 ! N #1 charge is located at N #1 ql(3,1) = 0.16536 ! N #2 charge is located at N #2 eps(1,1) = 1.00000 ! N #1 diameter / N diameter sig(1,1) = 1.00000 ! N #1 epsilon / N epsilon eps(2,1) = 1.00000 ! N #2 diameter / N diameter sig(2,1) = 1.00000 ! N #2 epsilon / N epsilon Na(1) = 2 ! # of LJ sites Nb(1,1) = num_mol(1)*num_mol(4)*num_mol(4) ! number of molecules of N2 Nb(1,2) = 0 ! # of N2 in phase (2) Nc(1) = 3 ! # of charge sites ELSE IF (N2_model.EQ.2) THEN ! This is a new N2 model ONLY to test against literature isotherms bl(1,1) = 0.16476 ! position of N #1 bl(2,1) = 0.16476 ! position of N #2 ql(1,1) = 0.000 ! charge #1 is located at COM of N2 ql(2,1) = 0.000 ! N #1 charge is located at N #1 ql(3,1) = 0.000 ! N #2 charge is located at N #2 eps(1,1) = 1.03846 ! N #1 diameter / N diameter sig(1,1) = 0.99940 ! N #1 epsilon / N epsilon eps(2,1) = 1.03846 ! N #2 diameter / N diameter sig(2,1) = 0.99940 ! N #2 epsilon / N epsilon Na(1) = 2 ! # of LJ sites Nb(1,1) = num_mol(1)*num_mol(4)*num_mol(4) ! number of molecules of N2 Nb(1,2) = 0 ! # of N2 in phase (2) Nc(1) = 0 ! # of charge sites ELSE IF (N2_model.EQ.3) THEN ! This is a TEST for chem pot. Nb(1,1) = num_mol(1)*num_mol(4)*num_mol(4) Nb(1,2) = 0 eps(1,1) = 1.000 ! model #1 (primary model) sig(1,1) = 0.769 Na(1) = 1 Nc(1) = 0 bl(1,1) = 0.0 ELSE; STOP END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! HYDROGEN ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Nb(2,1) = num_mol(2)*num_mol(4)*num_mol(4) ! number of molecules of N2 Nb(2,2) = 0 ! # of H2 in phase (2) ! eps(1,2) = 1.044 ! model #1 (primary model) ! sig(1,2) = 0.878 ! Na(2) = 1 ! Nc(2) = 0 ! bl(1,2) = 0.0 eps(1,2) = 1.008 ! model #2 (secondary model) sig(1,2) = 0.891 Na(2) = 2 Nc(2) = 3 bl(1,2) = 0.1116 bl(2,2) = 0.1116 ql(1,2) = 0.1116 ql(2,2) = 0.1116 ql(3,2) = 0.0000 charge(1,2) = 0.615 charge(2,2) = 0.615 charge(3,2) = -1.230 ! eps(1,2) = 0.9148 ! model #3 ! sig(1,2) = 0.8940 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! AMMONIA ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NH3_model = 2 IF (NH3_model.EQ.1) THEN ! 1-LJ model charge(1,3) = -1.026 ! partial charge on the nitrogen charge(2,3) = 0.342 ! partial charge on H #1 charge(3,3) = 0.342 ! partial charge on H #2 charge(4,3) = 0.342 ! partial charge on H #3 eps(1,3) = 2.9052 ! NH3 epsilon / N epsilon sig(1,3) = 1.0120 ! NH3 diameter / N diameter eps(2,3) = 0.0000 sig(2,3) = 1.0000 bl(1,3) = 0.000 ! nitrogen is located at COM of ammonia bl(2,3) = 1.000 ! dummy position to calculate orientation ql(1,3) = 0.000 ! nitrogen charge is located at COM of NH3 ql(2,3) = 0.30494 ! H#1 charge is located at ql/sigma(nitrogen=3.36ang) ql(3,3) = 0.30494 ! H#2 charge is located at ql/sigma(nitrogen=3.36ang) ql(4,3) = 0.30494 ! H#3 charge is located at ql/sigma(nitrogen=3.36ang) phi = 0.0 ! initial rotation about x-axis theta = 0.37696*pi ! angle made with x-axis, so that HNH angle=106.67 Na(3) = 2 ! # of LJ sites on each molecule Nb(3,1) = num_mol(3)*num_mol(4)*num_mol(4) ! # of molecules of Ammonia Nb(3,2) = 0 ! # of NH3 in phase (2) Nc(3) = 4 ! # of charges on each molecule ELSE IF (NH3_model.EQ.2) THEN ! 4-LJ model charge(1,3) = -1.455 ! partial charge on the nitrogen charge(2,3) = 0.485 ! partial charge on H #1 charge(3,3) = 0.485 ! partial charge on H #2 charge(4,3) = 0.485 ! partial charge on H #3 eps(1,3) = 1.0000 ! N of NH3 epsilon / N epsilon eps(2,3) = 0.5797 ! H #1 of NH3 eps(3,3) = 0.5797 ! H #2 of NH3 eps(4,3) = 0.5797 ! H #3 of NH3 sig(1,3) = 1.0000 ! N of NH3 diameter / N diameter sig(2,3) = 0.6777 sig(3,3) = 0.6777 sig(4,3) = 0.6777 bl(1,3) = 0.00000 ! N LJ site is located at COM of NH3 bl(2,3) = 0.30494 ! dummy position to calculate orientation bl(3,3) = 0.30494 bl(4,3) = 0.30494 ql(1,3) = 0.04699 ! nitrogen charge is offset towards the H's ql(2,3) = 0.30494 ! H#1 charge is located at ql/sigma(nitrogen=3.36ang) ql(3,3) = 0.30494 ! H#2 charge is located at ql/sigma(nitrogen=3.36ang) ql(4,3) = 0.30494 ! H#3 charge is located at ql/sigma(nitrogen=3.36ang) phi = 0.0 ! initial rotation about x-axis theta = 0.37696*pi ! angle made with x-axis Na(3) = 4 ! # of LJ sites on each molecule Nb(3,1) = num_mol(3)*num_mol(4)*num_mol(4) ! # of molecules of Ammonia Nb(3,2) = 0 ! # of NH3 in phase (2) Nc(3) = 4 ! # of charges on each molecule END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! NEW AMMONIA: 12.2.99 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This new model for ammonia places LJ sites on the ! ! hydrogens, so that I get a more realistic potential ! ! between the hydrogen and the carbon surface. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! charge(1,3) = -1.455 ! partial charge on the nitrogen ! charge(2,3) = 0.485 ! partial charge on H #1 (units of e charge) ! charge(3,3) = 0.485 ! partial charge on H #2 ! charge(4,3) = 0.485 ! partial charge on H #3 ! eps(1,3) = 1.000 ! N epsilon / N epsilon ! eps(2,3) = 0.5797 ! H #1 epsilon / N epsilon ! eps(3,3) = 0.5797 ! H #2 epsilon / N epsilon ! eps(4,3) = 0.5797 ! H #3 epsilon / N epsilon ! sig(1,3) = 1.000 ! N diameter / N diameter ! sig(2,3) = 0.6777 ! H #1 diameter / N diameter ! sig(3,3) = 0.6777 ! H #2 diameter / N diameter ! sig(4,3) = 0.6777 ! H #3 diameter / N diameter ! bl(1,3) = 0.000 ! nitrogen LJ site is located at COM of ammonia ! bl(2,3) = 0.30494 ! H#1 LJ site is located at COM of ammonia ! bl(3,3) = 0.30494 ! H#2 LJ site is located at COM of ammonia ! bl(4,3) = 0.30494 ! H#3 LJ site is located at COM of ammonia ! ql(1,3) = 0.04699 ! nitrogen charge site is offset from COM of NH3 ! ql(2,3) = 0.30494 ! H#1 charge is located at ql/sigma(nitrogen=3.32ang) ! ql(3,3) = 0.30494 ! H#2 charge is located at ql/sigma(nitrogen=3.32ang) ! ql(4,3) = 0.30494 ! H#3 charge is located at ql/sigma(nitrogen=3.32ang) ! phi = 0.0 ! initial rotation about x-axis ! theta = 0.37696*pi ! angle made with x-axis, so that HNH angle=106.67 ! Na(3) = 4 ! # of LJ sites on each molecule ! Nb(3,1) = num_mol(3)*num_mol(4)*num_mol(4) ! # of molecules of Ammonia ! Nb(3,2) = 0 ! # of NH3 in phase (2) ! Nc(3) = 4 ! # of charges on each molecule !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! METHANE ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(inert.EQ.1) THEN IF(kk.EQ.1) Nb(4,1) = 91 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(4,1) = 84 IF(kk.EQ.3) Nb(4,1) = 81 IF(kk.EQ.4) Nb(4,1) = 79 IF(kk.EQ.5) Nb(4,1) = 78 IF(kk.EQ.6) Nb(4,1) = 76 IF(kk.EQ.7) Nb(4,1) = 75 IF(kk.EQ.8) Nb(4,1) = 74 IF(kk.EQ.9) Nb(4,1) = 74 IF(kk.EQ.10) Nb(4,1) = 74 END IF IF(inert.EQ.2) THEN IF(kk.EQ.1) Nb(4,1) = 108 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(4,1) = 99 IF(kk.EQ.3) Nb(4,1) = 93 IF(kk.EQ.4) Nb(4,1) = 89 IF(kk.EQ.5) Nb(4,1) = 87 IF(kk.EQ.6) Nb(4,1) = 84 IF(kk.EQ.7) Nb(4,1) = 82 IF(kk.EQ.8) Nb(4,1) = 81 IF(kk.EQ.9) Nb(4,1) = 80 IF(kk.EQ.10) Nb(4,1) = 79 END IF IF(inert.EQ.3) THEN IF(kk.EQ.1) Nb(4,1) = 108 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(4,1) = 101 IF(kk.EQ.3) Nb(4,1) = 96 IF(kk.EQ.4) Nb(4,1) = 93 IF(kk.EQ.5) Nb(4,1) = 89 IF(kk.EQ.6) Nb(4,1) = 87 IF(kk.EQ.7) Nb(4,1) = 84 IF(kk.EQ.8) Nb(4,1) = 82 IF(kk.EQ.9) Nb(4,1) = 81 IF(kk.EQ.10) Nb(4,1) = 79 END IF IF(inert.EQ.4) THEN IF(kk.EQ.1) Nb(4,1) = 115 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(4,1) = 111 IF(kk.EQ.3) Nb(4,1) = 108 IF(kk.EQ.4) Nb(4,1) = 105 IF(kk.EQ.5) Nb(4,1) = 102 IF(kk.EQ.6) Nb(4,1) = 100 IF(kk.EQ.7) Nb(4,1) = 97 IF(kk.EQ.8) Nb(4,1) = 95 IF(kk.EQ.9) Nb(4,1) = 93 IF(kk.EQ.10) Nb(4,1) = 91 END IF Nb(4,2) = 0 Na(4) = 1 Nc(4) = 0 eps(1,4) = 4.0714 sig(1,4) = 1.1497 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ARGON ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(inert.EQ.1) THEN IF(kk.EQ.1) Nb(5,1) = 30 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(5,1) = 28 IF(kk.EQ.3) Nb(5,1) = 27 IF(kk.EQ.4) Nb(5,1) = 26 IF(kk.EQ.5) Nb(5,1) = 26 IF(kk.EQ.6) Nb(5,1) = 25 IF(kk.EQ.7) Nb(5,1) = 25 IF(kk.EQ.8) Nb(5,1) = 25 IF(kk.EQ.9) Nb(5,1) = 25 IF(kk.EQ.10) Nb(5,1) = 25 END IF IF(inert.EQ.2) THEN IF(kk.EQ.1) Nb(5,1) = 36 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(5,1) = 33 IF(kk.EQ.3) Nb(5,1) = 31 IF(kk.EQ.4) Nb(5,1) = 30 IF(kk.EQ.5) Nb(5,1) = 29 IF(kk.EQ.6) Nb(5,1) = 28 IF(kk.EQ.7) Nb(5,1) = 27 IF(kk.EQ.8) Nb(5,1) = 27 IF(kk.EQ.9) Nb(5,1) = 27 IF(kk.EQ.10) Nb(5,1) = 26 END IF IF(inert.EQ.3) THEN IF(kk.EQ.1) Nb(5,1) = 36 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(5,1) = 34 IF(kk.EQ.3) Nb(5,1) = 32 IF(kk.EQ.4) Nb(5,1) = 31 IF(kk.EQ.5) Nb(5,1) = 30 IF(kk.EQ.6) Nb(5,1) = 29 IF(kk.EQ.7) Nb(5,1) = 28 IF(kk.EQ.8) Nb(5,1) = 27 IF(kk.EQ.9) Nb(5,1) = 27 IF(kk.EQ.10) Nb(5,1) = 26 END IF IF(inert.EQ.4) THEN IF(kk.EQ.1) Nb(5,1) = 38 ! on basis of 1728 molecules total IF(kk.EQ.2) Nb(5,1) = 37 IF(kk.EQ.3) Nb(5,1) = 36 IF(kk.EQ.4) Nb(5,1) = 35 IF(kk.EQ.5) Nb(5,1) = 34 IF(kk.EQ.6) Nb(5,1) = 33 IF(kk.EQ.7) Nb(5,1) = 32 IF(kk.EQ.8) Nb(5,1) = 32 IF(kk.EQ.9) Nb(5,1) = 31 IF(kk.EQ.10) Nb(5,1) = 30 END IF Nb(5,2) = 0 Na(5) = 1 Nc(5) = 0 eps(1,5) = 3.2912 sig(1,5) = 1.0256 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Steele Potential Parameters ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sig_c = 1.024 eps_c = 0.769231 DO i=1,nspecies DO j=1,Na(i) st_sig(j,i) = 0.5*(sig_c+sig(j,i)) ! sigma for Steele Potential (LB mixing) st_eps(j,i) = SQRT(eps_c*eps(j,i)) ! epsilon for Steele Potential (LB mixing) END DO END DO delta=1.009 ! graphite layer spacing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Lorentz-Berthelot Mixing Rules ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO i=1,nspecies DO j=1,nspecies atom1: DO k=1,Na(i) atom2: DO l=1,Na(j) eps_m(k,i,l,j) = SQRT(eps(k,i) * eps(l,j)) ! epsilon b/t type(i) and (j) sig_m(k,i,l,j) = 0.50*(sig(k,i) + sig(l,j)) ! sigma b/t type(i) and (j) eps_m(l,j,k,i) = eps_m(l,j,k,i) ! should be no ordering effect sig_m(l,j,k,i) = sig_m(l,j,k,i) ! should be no ordering effect END DO atom2 END DO atom1 END DO END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! r=0.0;xr=0.0;yr=0.0;zr=0.0 t=0.0;loop=0 numave=0 nsub=1 cumfracA1=0.0 cumfracA2=0.0 cumfracT=0.0 cumNb_b=0.0 cumNb_s=0.0 cumN2_b=0.0 cumN2_s=0.0 cumH2_b=0.0 cumH2_s=0.0 cumU_b=0.0 cumU_s=0.0 cum_P=0.0 cum_vol=0.0 Nstatus=0 chem_pot=0.0 cum_chem_pot=0.0 ! histC1 = 0.0 ! RDF data arrays ! histC2 = 0.0 ndel=0; nins=0 ! counters for # of particle swaps b/t phases acc1=0.0; acc2=0.0; acc3_b=0.0; acc3_s=0.0; acc4_b=0.0; acc4_s=0.0 att1=0.0; att2=0.0; att3_b=0.0; att3_s=0.0; att4_b=0.0; att4_s=0.0 sd1=0.0;sd2=0.0;sd3=0.0;sd4=0.0;sd5=0.0;sd6=0.0;sd7=0.0;sd8=0.0 rho_m(1) = Rho_set(kk) ! initialize density REWIND(100) ! rewind the data files for next T/P conditions REWIND(110) REWIND(400) REWIND(1400) IF (nanim.EQ.1) REWIND(900) IF (nanim.EQ.1) REWIND(1000) IF((Niso.EQ.1).AND.(kk.GT.TP_i)) GOTO 500 Vol(1) = REAL(num_mol(5))/rho_m(1) ! Volume of bulk Vol0 = Vol(1) Vol(2) = Vol(1)*V_pore/1000.0 ! Volume of pore box(1,1) = Vol(1)**(1.0/3.0) box(2,1) = Vol(1)**(1.0/3.0) box(3,1) = Vol(1)**(1.0/3.0) box(1,2) = SQRT(Vol(2)/box(3,2)) box(2,2) = SQRT(Vol(2)/box(3,2)) DO k=1,num_mol(1),1; DO j=1,num_mol(4),1; DO i=1,num_mol(4),1 ! Generate initial N2 positions xm(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),1,1) = (REAL(i)-0.5) ym(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),1,1) = (REAL(j)-0.5) zm(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),1,1) = (REAL(k)-0.5) END DO; END DO ; END DO DO k=1,num_mol(2),1; DO j=1,num_mol(4),1; DO i=1,num_mol(4),1 ! Generate initial H2 positions xm(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),2,1) = (REAL(i)-0.5) ym(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),2,1) = (REAL(j)-0.5) zm(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),2,1) = (REAL(k+num_mol(1))-0.5) END DO; END DO ; END DO DO k=1,num_mol(3),1; DO j=1,num_mol(4),1; DO i=1,num_mol(4),1 ! Generate initial NH3 positions xm(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),3,1) = (REAL(i)-0.5) ym(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),3,1) = (REAL(j)-0.5) zm(i+num_mol(4)*(j-1)+num_mol(4)*num_mol(4)*(k-1),3,1) = (REAL(k+num_mol(1)+num_mol(2))-0.5) END DO; END DO ; END DO IF(cut.GT.(box(1,1)/2.0)) THEN ! check the length of the cutoff PRINT *,'LJ Cut-off radius is too high!' PRINT *,'LJ Cut-off: ',cut PRINT *,'Box length/2: ',box(1,1)/2.0 READ * END IF IF(el_cut.GT.(box(1,1)/2.0)) THEN ! check the length of the cutoff PRINT *,'Electrostatic Cut-off radius is too high!' PRINT *,'Electrostatic Cut-off: ',el_cut PRINT *,'Box length/2: ',box(1,1)/2.0 READ * END IF DO k=1,nspecies DO i=1,Nb(k,1) ! Scale to initial density xm(i,k,1) = xm(i,k,1)*box(1,1)/REAL(num_mol(4)) ym(i,k,1) = ym(i,k,1)*box(2,1)/REAL(num_mol(4)) zm(i,k,1) = zm(i,k,1)*box(3,1)/REAL(num_mol(4)) xm(i,k,1) = xm(i,k,1) - (box(1,1)/2.0) ym(i,k,1) = ym(i,k,1) - (box(2,1)/2.0) zm(i,k,1) = zm(i,k,1) - (box(3,1)/2.0) END DO END DO DO i=1,Nb(1,1) ! Initial orientation for N2 ex(1,i,1,1) = 1.0 ey(1,i,1,1) = 0.0 ez(1,i,1,1) = 0.0 ex(2,i,1,1) = -1.0 ey(2,i,1,1) = 0.0 ez(2,i,1,1) = 0.0 qx(1,i,1,1) = 0.0 qy(1,i,1,1) = 0.0 qz(1,i,1,1) = 0.0 qx(2,i,1,1) = 1.0 qy(2,i,1,1) = 0.0 qz(2,i,1,1) = 0.0 qx(3,i,1,1) = -1.0 qy(3,i,1,1) = 0.0 qz(3,i,1,1) = 0.0 END DO DO i=1,Nb(2,1) ! Initial orientation for H2 ex(1,i,1,1) = 1.0 ey(1,i,1,1) = 0.0 ez(1,i,1,1) = 0.0 ex(2,i,1,1) = -1.0 ey(2,i,1,1) = 0.0 ez(2,i,1,1) = 0.0 qx(1,i,1,1) = 0.0 qy(1,i,1,1) = 0.0 qz(1,i,1,1) = 0.0 qx(2,i,1,1) = 1.0 qy(2,i,1,1) = 0.0 qz(2,i,1,1) = 0.0 qx(3,i,1,1) = -1.0 qy(3,i,1,1) = 0.0 qz(3,i,1,1) = 0.0 END DO DO i=1,Nb(2,1) ! Initial orientation for H2 ex(1,i,2,1) = 1.0 ey(1,i,2,1) = 0.0 ez(1,i,2,1) = 0.0 ex(2,i,2,1) = -1.0 ey(2,i,2,1) = 0.0 ez(2,i,2,1) = 0.0 qx(1,i,2,1) = 1.0 qy(1,i,2,1) = 0.0 qz(1,i,2,1) = 0.0 qx(2,i,2,1) = -1.0 qy(2,i,2,1) = 0.0 qz(2,i,2,1) = 0.0 qx(3,i,2,1) = 0.0 qy(3,i,2,1) = 0.0 qz(3,i,2,1) = 0.0 END DO DO i=1,Nb(3,1) ! Initial orientation for NH3 qx(1,i,3,1) = 1.0 ! Assign initial charge orientation qx(2,i,3,1) = COS(theta) qx(3,i,3,1) = COS(theta) qx(4,i,3,1) = COS(theta) qy(1,i,3,1) = 0.0 qy(2,i,3,1) = SIN(theta)*COS(0.0) qy(3,i,3,1) = SIN(theta)*COS(0.667*pi) qy(4,i,3,1) = SIN(theta)*COS(-0.667*pi) qz(1,i,3,1) = 0.0 qz(2,i,3,1) = SIN(theta)*SIN(0.0) qz(3,i,3,1) = SIN(theta)*SIN(0.667*pi) qz(4,i,3,1) = SIN(theta)*SIN(-0.667*pi) ex(1,i,3,1) = 1.0 ! Assign initial LJ orientation ex(2,i,3,1) = COS(theta) ex(3,i,3,1) = COS(theta) ex(4,i,3,1) = COS(theta) ey(1,i,3,1) = 0.0 ey(2,i,3,1) = SIN(theta)*COS(0.0) ey(3,i,3,1) = SIN(theta)*COS(0.667*pi) ey(4,i,3,1) = SIN(theta)*COS(-0.667*pi) ez(1,i,3,1) = 0.0 ez(2,i,3,1) = SIN(theta)*SIN(0.0) ez(3,i,3,1) = SIN(theta)*SIN(0.667*pi) ez(4,i,3,1) = SIN(theta)*SIN(-0.667*pi) END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! INSERT INERT MOLECULES ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Some of the Nitrogen molecules will be replaced by Argon, and ! ! some of the Hydrogen molecules will be replaced by methane ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO i=Nb(1,1)-Nb(5,1)+1, Nb(1,1) ! replace N2 with Ar xm((i-Nb(1,1)+Nb(5,1)),5,1) = xm(i,1,1) ym((i-Nb(1,1)+Nb(5,1)),5,1) = ym(i,1,1) zm((i-Nb(1,1)+Nb(5,1)),5,1) = zm(i,1,1) END DO Nb(1,1) = Nb(1,1) - Nb(5,1) DO i=Nb(2,1)-Nb(4,1)+1, Nb(2,1) ! replace H2 with CH4 xm((i-Nb(2,1)+Nb(4,1)),4,1) = xm(i,2,1) ym((i-Nb(2,1)+Nb(4,1)),4,1) = ym(i,2,1) zm((i-Nb(2,1)+Nb(4,1)),4,1) = zm(i,2,1) END DO Nb(2,1) = Nb(2,1) - Nb(4,1) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Uold(1) = Force_calc(xm,ym,zm,ex,ey,ez,qx,qy,qz,1) ! Initial values for Vir and PE Uold(2) = 0.0 Vir = VirNew Press = P_temp Vol(2) = V_pore*1000.0 ! Define the volume of the pore phase box(1,2) = SQRT(Vol(2)/box(3,2)) box(2,2) = SQRT(Vol(2)/box(3,2)) 500 CONTINUE ! Skip to here for iso. calc. PRINT 230,'Temperature: ',Temp(kk)*36.4 ! Show info to screen PRINT 230,'Pressure: ',P_set(kk)*137.334 PRINT 230,'Molecular density: ',Rho_set(kk) ! PRINT *,'Pore width: ',box(3,2) PRINT *,'Run #: ',num_run WRITE (300,230) 'Temperature/K: ',Temp(kk)*36.4 ! Save initial info WRITE (300,230) 'Pressure/bar: ',P_set(kk)*137.334 WRITE (300,230) 'Density/(N/sigma^3): ',Rho_set(kk) ! WRITE (300,*) 'Pore height: ',box(3,2) 200 FORMAT (1x,I4) 210 FORMAT (1x,A10) 220 FORMAT (1x,A3,3F9.4) 230 FORMAT (1x,A,F8.2) END SUBROUTINE Initialize !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Show_status INCLUDE 'ammonia.inc' DIMENSION :: Uold2(N) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! calc1(s): average energy per molecule ! ! calc2(s): mole fraction of monomer ! ! calc3(s): total number of atoms ! ! ncalc3(s): molecular density ! ! calc4(s): atomic density ! ! ! ! ratio4: reaction acceptance ratio ! ! ratio5: unreaction acceptance ratio ! ! ratio6: move acceptance ratio ! ! ratio7: rotation acceptance ratio ! ! ratio8: volume change accpetance ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! nTot_mol_1 = 0 nTot_mol_2 = 0 ratio8 = 0.0 IF(Nrdf.EQ.1) THEN Uold2(1) = Force_calc(xm,ym,zm,ex,ey,ez,qx,qy,qz,1) Press = P_temp Vir = VirNew END IF Uold2(2) = Force_calc(xm,ym,zm,ex,ey,ez,qx,qy,qz,2) count: DO i=1,nspecies ! Count # of mol in each phase nTot_mol_1 = nTot_mol_1 + Nb(i,1) nTot_mol_2 = nTot_mol_2 + Nb(i,2) END DO count IF(nTot_mol_1.GT.0) THEN ! Phase #1 system proprties calc1 = Uold(1)/REAL(nTot_mol_1) calc2 = REAL(Nb(3,1))/REAL(nTot_mol_1) ! ncalc3 = (2*Nb(1)+Na(1)) calc4 = REAL(nTot_mol_1)/Vol(1) ! calc5 = REAL(Na(1)+2*Nb(1))/Vol(1) END IF ratio4=0.0;ratio5=0.0;ratio8=0.0 ratio6_b=0.0; ratio6_s=0.0 ratio7_b=0.0; ratio7_s=0.0 IF (att1.GT.1.0) ratio4 = acc1/att1 IF (att2.GT.1.0) ratio5 = acc2/att2 IF (att3_b.GT.1.0) ratio6_b = acc3_b/att3_b IF (att3_s.GT.1.0) ratio6_s = acc3_s/att3_s IF (att4_b.GT.1.0) ratio7_b = acc4_b/att4_b IF (att4_s.GT.1.0) ratio7_s = acc4_s/att4_s IF (att5.GT.1.0) ratio8 = acc5/att5 IF(nTot_mol_2.GT.0) THEN ! Phase #2 system properties calc1s = Uold(2)/REAL(nTot_mol_2) calc2s = REAL(Nb(3,2))/REAL(nTot_mol_2) calc4s = REAL(nTot_mol_2)/(Vol(2)*(box(3,2)-sig_c)/box(3,2)) END IF calc2t = REAL(Nb(3,2)+Nb(3,1))/REAL(nTot_mol_1 + nTot_mol_2) IF(nprint.EQ.1) THEN ! nprint = 1 for screen printout PRINT *,'________________________________' PRINT *,'Bulk energy (full)',Uold2(1) PRINT *,'Bulk energy (part)',Uold(1) PRINT *,'Difference: ',(Uold2(1) - Uold(1)) PRINT *,'Slit energy (full)',Uold2(2) PRINT *,'Slit energy (part)',Uold(2) PRINT *,'Difference: ',Uold2(2) - Uold(2) PRINT * PRINT 230,'Time: ',t PRINT 100,'Temperature(K): ',Temp(kk)*36.4 ! convert to real units PRINT 100,'Temperature(C): ',Temp(kk)*36.4-273.15 ! convert to real units PRINT 100,'Pore width (sigma): ',box(3,2) PRINT 100,'Energy (bulk): ',calc1 PRINT 100,'Energy (slit): ',calc1s PRINT 101,'Mols in bulk (1,2,3,4,5): ',Nb(1,1),Nb(2,1),Nb(3,1),Nb(4,1),Nb(5,1) PRINT 101,'Mols in pore (1,2,3,4,5): ',Nb(1,2),Nb(2,2),Nb(3,2),Nb(4,2),Nb(5,2) PRINT 302,'Total # of atoms: ',2*Nb(1,1)+2*Nb(2,1)+4*Nb(3,1)+2*Nb(1,2)+2*Nb(2,2)+4*Nb(3,2) PRINT 100,'Mole frac ammonia (bulk): ',calc2 PRINT 100,'Mole frac ammonia (slit): ',calc2s PRINT *,'Total# of mols (bulk): ',nTot_mol_1 PRINT *,'Total# of mols (slit): ',nTot_mol_2 PRINT 240,'Volume (bulk): ',Vol(1) PRINT 100,'Volume (slit): ',Vol(2) PRINT 300,'Pressure set (bars): ',P_set(kk)*137.334 ! convert to real units PRINT 300,'Pressure (bars): ',Press*137.334 PRINT 100,'Molecular density (bulk): ',calc4 PRINT 100,'Molecular density (slit): ',calc4s PRINT * PRINT 100,'Rot. acceptance ratio (b): ',ratio7_b PRINT 100,'Rot. acceptance ratio (s): ',ratio7_s PRINT 100,'Move acceptance ratio (b): ',ratio6_b PRINT 100,'Move acceptance ratio (s): ',ratio6_s PRINT 100,'Maximum displacement (bulk):',disp_b PRINT 100,'Maximum displacement (slit):',disp_s PRINT 100,'Maximum rotation (bulk): ',(Dgamax_b*180.0/pi) PRINT 100,'Maximum rotation (slit): ',(Dgamax_s*180.0/pi) PRINT 100,'Reaction acceptance ratio: ',ratio4 PRINT 100,'Unreaction acceptance ratio:',ratio5 PRINT 100,'Volume change acceptance: ',ratio8 PRINT 100,'Maximum volume change ',dbox PRINT *,'________________________________' END IF !!! Write energy, pressure, Xa, and volume results to a file !!! WRITE (110,260) t,Press*137.334,calc1 ! write: pressure, E(bulk) WRITE (400,250) t,calc2 ! write: x_NH3(bulk) ! WRITE (1400,*) t,Nb(1,2),Nb(2,2),Nb(3,2) ! write: #N2, #H2, #NH3 (pore) !!! Make adjustments in move parameters disp_b = disp_b/(1.4-ratio6_b) ! adjust max displacement in bulk IF (disp_b.GT.(box(1,1)/10.0)) disp_b = box(1,1)/10.0 IF (disp_b.LT.0.10) disp_b = 0.10 disp_s = disp_s/(1.4-ratio6_s) ! adjust max displacement in pore IF (disp_s.GT.(box(3,2)/2.0)) disp_s = box(3,2)/2.0 IF (disp_s.LT.0.1) disp_s = 0.1 Dgamax_b = Dgamax_b/(1.4 - ratio7_b) ! max rotation in bulk IF (Dgamax_b.GT.pi) Dgamax_b = pi IF (Dgamax_b.LT.(pi/8.0)) Dgamax_b = (pi/8.0) Dgamax_s = Dgamax_s/(1.4 - ratio7_s) ! max rotation in pore IF (Dgamax_s.GT.pi) Dgamax_s = pi IF (Dgamax_s.LT.(pi/8.0)) Dgamax_s = (pi/8.0) dbox = dbox/(1.40 - ratio8) ! adjust max change in volume IF (dbox.GT.0.2*box(1,1)) dbox = 0.25*box(1,1) IF (dbox.LT.1.0e-03*box(1,1)) dbox = 1.0e-03*box(1,1) ! acc1 = 0.0 ! zero counters ! att1 = 0.0 ! acc2 = 0.0 ! att2 = 0.0 acc3_b = 0.0 acc3_s = 0.0 att3_b = 0.0 att3_s = 0.0 acc4_b = 0.0 acc4_s = 0.0 att4_b = 0.0 att4_s = 0.0 Nstatus = Nstatus + 1 ! count # of loops through 'Status' subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 100 FORMAT (1x,A,F12.5) 101 FORMAT (1x,A,5(1x,I4)) 200 FORMAT (1x,I4) 210 FORMAT (1x,A10) 220 FORMAT (1x,A3,3F9.4) 230 FORMAT (1x,A,F12.0) 240 FORMAT (1x,A,F12.3) 250 FORMAT (1x,F12.2,2F10.5) 260 FORMAT (1x,F12.2,2x,F12.4,3F12.4) 300 FORMAT (1x,A,F12.6) 301 FORMAT (1x,A,3I5) 302 FORMAT (1x,A,I5) END SUBROUTINE Show_status !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Move !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Thus subroutine will displace a randomly selected molecule in any ! ! phase(nph). The COM will be moved in this subroutine and then the ! ! position and orientation vectors will be passed to BForce, which ! ! will calculate the atomic positions and calculate the change in ! ! potential E for the move. ! ! ! ! Tot_mol_nph: total m# of molecules in phase (nph) ! ! ntype: the species type of molecule to move ! ! numpart: the molecule # in the array to move ! ! orient1: orientations of atoms in displaced molecule ! ! orient2: orientations of charges in displaced molecule ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION :: orient1(natoms,3),orient2(ncgs,3),coord(3) movea = 0 moveb = 0 deltaU = 0.0 NU1accept = 0 dpos = 0.0 Tot_mol_nph = 0 orient = 0.0 coord = 0.0 dVir = 0.0 DO i=1,nspecies Tot_mol_nph = Tot_mol_nph + Nb(i,nph) END DO IF(Tot_mol_nph.LE.0) RETURN ! check for mols in phase (nph) SELECT CASE (nph) ! find max. disp. in phase (nph) CASE(1) dpos = disp_b att3_b = att3_b + 1.0 CASE(2) dpos = disp_s att3_s = att3_s + 1.0 END SELECT numpart = INT(AINT(ran3(iseed)*REAL(Tot_mol_nph))) + 1 ! pick a random particle DO i=1,nspecies ! pick the exact mol to move ntype = i IF (numpart.LE.Nb(i,nph)) EXIT numpart=numpart-Nb(i,nph) END DO coord(1) = xm(numpart,ntype,nph) + dpos*(ran3(iseed) - 0.5) ! new x-position coord(2) = ym(numpart,ntype,nph) + dpos*(ran3(iseed) - 0.5) ! new y-position coord(3) = zm(numpart,ntype,nph) + dpos*(ran3(iseed) - 0.5) ! new z-position DO i=1,Na(ntype) ! save atom orientations orient1(i,1) = ex(i,numpart,ntype,nph) orient1(i,2) = ey(i,numpart,ntype,nph) orient1(i,3) = ez(i,numpart,ntype,nph) END DO DO i=1,Nc(ntype) ! save charge orientations orient2(i,1) = qx(i,numpart,ntype,nph) orient2(i,2) = qy(i,numpart,ntype,nph) orient2(i,3) = qz(i,numpart,ntype,nph) END DO deltaU = BForce(orient1,orient2,coord,numpart,ntype,nph) ! calc. change in U IF (novr.EQ.1) RETURN ! return to main loop if an overlap is created IF ((deltaU/Temp(kk)).GT.70.0) RETURN ! IF (deltaU.LE.0.0) THEN ! acceptance probability NU1accept = 1 ELSE z1 = MIN(1.0,EXP(-deltaU/Temp(kk))) z2 = ran3(iseed)*1.0 IF (z2.LT.z1) NU1accept = 1 END IF acceptance: IF (NU1accept.EQ.1) THEN ! update position and U SELECT CASE (nph) CASE (1); acc3_b = acc3_b + 1.0 CASE (2); acc3_s = acc3_s + 1.0 END SELECT Uold(nph) = Uold(nph) + deltaU xm(numpart,ntype,nph) = coord(1) ym(numpart,ntype,nph) = coord(2) zm(numpart,ntype,nph) = coord(3) END IF acceptance END SUBROUTINE Move !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Rotate !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This performs rotations in either phase (nph). Only the center of ! ! mass and orientation of the rotated molecule is passed to the force ! ! subroutine (BForce). If accepted, the orientation is updated at ! ! the end. However, the atomic positions are not explicitly updated. ! ! The center of mass of the molecule is orientation vector '1'. ! ! Therefore, this index is skipped since rotations will have no ! ! effect on the center of mass. If this index IS included, it will ! ! cause an overflow when it is normalized in the rotation loop. ! ! This should also save time. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION :: orient1(natoms,3),orient2(ncgs,3),coord(3) NU2accept = 0 rnorm = 0.0 Drot = 0.0 orient = 0.0 coord = 0.0 nTot_mol_nph = 0 deltaU = 0.0 DO i=1,nspecies nTot_mol_nph = nTot_mol_nph + Nb(i,nph) END DO IF(nTot_mol_nph.LE.0) RETURN ! check for mols in phase (nph) SELECT CASE (nph) ! Choose the phase (nph=1: bulk phase, nph=2: pore) CASE(1) Drot = Dgamax_b att4_b = att4_b + 1.0 CASE(2) Drot = Dgamax_s att4_s = att4_s + 1.0 END SELECT numpart = INT(AINT(ran3(iseed)*REAL(nTot_mol_nph))) + 1 DO i=1,nspecies ! pick the exact mol to rotate ntype = i IF (numpart.LE.Nb(i,nph)) EXIT numpart=numpart-Nb(i,nph) END DO IF(ntype.EQ.2) RETURN ! return if a H2 is selected (1-LJ has no orientation) Iaxis= INT(ran3(iseed)*3.0) + 1 ! pick an axis Dgamma= (ran3(iseed) - 0.5) * Drot ! random orientation cosDg = COS(Dgamma) ! set up matrix sinDg = SIN(Dgamma) SELECT CASE (Iaxis) ! choose the axis of rotation CASE(1) ! rotate about x-axis DO i=1,Na(ntype) ! rotate atoms orient1(i,1) = ex(i,numpart,ntype,nph) orient1(i,2) = cosDg * ey(i,numpart,ntype,nph) + sinDg * ez(i,numpart,ntype,nph) orient1(i,3) = cosDg * ez(i,numpart,ntype,nph) - sinDg * ey(i,numpart,ntype,nph) END DO DO i=1,Nc(ntype) ! rotate charges orient2(i,1) = qx(i,numpart,ntype,nph) orient2(i,2) = cosDg * qy(i,numpart,ntype,nph) + sinDg * qz(i,numpart,ntype,nph) orient2(i,3) = cosDg * qz(i,numpart,ntype,nph) - sinDg * qy(i,numpart,ntype,nph) END DO CASE(2) ! rotate about y-axis DO i=1,Na(ntype) ! rotate atoms orient1(i,1) = cosDg * ex(i,numpart,ntype,nph) - sinDg * ez(i,numpart,ntype,nph) orient1(i,2) = ey(i,numpart,ntype,nph) orient1(i,3) = cosDg * ez(i,numpart,ntype,nph) + sinDg * ex(i,numpart,ntype,nph) END DO DO i=1,Nc(ntype) ! rotate charges orient2(i,1) = cosDg * qx(i,numpart,ntype,nph) - sinDg * qz(i,numpart,ntype,nph) orient2(i,2) = qy(i,numpart,ntype,nph) orient2(i,3) = cosDg * qz(i,numpart,ntype,nph) + sinDg * qx(i,numpart,ntype,nph) END DO CASE(3) ! rotate about z-axis DO i=1,Na(ntype) ! rotate atoms orient1(i,1) = cosDg * ex(i,numpart,ntype,nph) + sinDg * ey(i,numpart,ntype,nph) orient1(i,2) = cosDg * ey(i,numpart,ntype,nph) - sinDg * ex(i,numpart,ntype,nph) orient1(i,3) = ez(i,numpart,ntype,nph) END DO DO i=1,Nc(ntype) ! rotate charges orient2(i,1) = cosDg * qx(i,numpart,ntype,nph) + sinDg * qy(i,numpart,ntype,nph) orient2(i,2) = cosDg * qy(i,numpart,ntype,nph) - sinDg * qx(i,numpart,ntype,nph) orient2(i,3) = qz(i,numpart,ntype,nph) END DO END SELECT DO i=1,Na(ntype) ! normalize the atom position unit vector IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(i,1)**2 + orient1(i,2)**2 + orient1(i,3)**2) orient1(i,1) = orient1(i,1)*rnorm orient1(i,2) = orient1(i,2)*rnorm orient1(i,3) = orient1(i,3)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge position unit vector IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(i,1)**2 + orient2(i,2)**2 + orient2(i,3)**2) orient2(i,1) = orient2(i,1)*rnorm orient2(i,2) = orient2(i,2)*rnorm orient2(i,3) = orient2(i,3)*rnorm END DO coord(1) = xm(numpart,ntype,nph) ! x-COM, position vector for BForce coord(2) = ym(numpart,ntype,nph) ! y-COM coord(3) = zm(numpart,ntype,nph) ! z-COM deltaU = BForce(orient1,orient2,coord,numpart,ntype,nph) ! Calc. U IF (novr.EQ.1) RETURN ! return if overlap occurs IF (deltaU/Temp(kk).GT.70.0) RETURN ! return if huge U increase IF (deltaU.LE.0.0) THEN ! accept if U decreases NU2accept = 1 ELSE ! acceptance criteria z1 = MIN(1.0,EXP(-deltaU/Temp(kk))) z2 = ran3(iseed)*1.0 IF (z2.LT.z1) NU2accept = 1 END IF acceptance: IF (NU2accept.EQ.1) THEN ! update orientation and U SELECT CASE (nph) CASE (1); acc4_b = acc4_b + 1.0 CASE (2); acc4_s = acc4_s + 1.0 END SELECT Uold(nph) = Uold(nph) + deltaU DO i=1,Na(ntype) ! update atomic orientations IF(ABS(bl(i,ntype)).LE.0.02) CYCLE ex(i,numpart,ntype,nph) = orient1(i,1) ey(i,numpart,ntype,nph) = orient1(i,2) ez(i,numpart,ntype,nph) = orient1(i,3) END DO DO i=1,Nc(ntype) ! update charge orientations IF(ABS(ql(i,ntype)).LE.0.02) CYCLE qx(i,numpart,ntype,nph) = orient2(i,1) qy(i,numpart,ntype,nph) = orient2(i,2) qz(i,numpart,ntype,nph) = orient2(i,3) END DO END IF acceptance END SUBROUTINE Rotate !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Forward !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This subroutine will perform a forward reaction step. For the ammonia ! ! synthesis reaction, one molecule of H2 and one molecule of N2 will each ! ! be changed to a molecule of NH3. In addition, two molecules of H2 will ! ! be deleted. The net result is the reaction N2 + 3H2 => 2NH3. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! q_ins(x/y/z,point_chg#,insertion#): inserted molecular charge orientation ! ! e_ins(x/y/z,LJ-site#,insertion#) : inserted LJ orientation relative to com ! ! com_ins(x/y/z,insertion#) : com of inseerted molecule ! ! n_rxn_part(mol#) : identity of molecule deleted/changed ! ! n_rxn_type(mol#) : type of molecule deleted/changed ! ! nnrxn_type : type of molecule to be inserted ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION com_ins(3,2),q_ins(3,ncgs,2),e_ins(3,natoms,2),n_rxn_part(4),n_rxn_type(4) DIMENSION orient1(3,natoms,2),orient2(3,ncgs,2),nn_rxn_type(2),num_h2(4),Iaxis(10) NU2accept=0 IF (t.GT.t_equil) att1=att1+1.0 IF(Nb(1,nph).LT.1) RETURN ! make sure N2 mols are present for rxn IF(Nb(2,nph).LT.3) RETURN ! make sure H2 mols are present for rxn numpart_n2 = INT(AINT(ran3(iseed)*REAL(Nb(1,nph)))) + 1 ! pick an N2 at random 1 numpart_h2a = INT(AINT(ran3(iseed)*REAL(Nb(2,nph)))) + 1 ! pick an H2 at random 2 numpart_h2b = INT(AINT(ran3(iseed)*REAL(Nb(2,nph)))) + 1 ! pick a different H2 at random IF(numpart_h2b.EQ.numpart_h2a) GOTO 2 3 numpart_h2c = INT(AINT(ran3(iseed)*REAL(Nb(2,nph)))) + 1 ! pick a different H2 at random IF((numpart_h2c.EQ.numpart_h2b).OR.(numpart_h2c.EQ.numpart_h2a)) GOTO 3 nins_rxn = 2; ndel_rxn = 4 n_rxn_part(1)=numpart_n2 ! Specify identity of N2 deleted n_rxn_part(2)=numpart_h2a ! Specify identity of H2 deleted n_rxn_part(3)=numpart_h2b ! Specify identity of H2 deleted n_rxn_part(4)=numpart_h2c ! Specify identity of H2 deleted n_rxn_type(1) = 1 ! Specify as N2 n_rxn_type(2) = 2 ! Specify as H2 n_rxn_type(3) = 2 ! Specify as H2 n_rxn_type(4) = 2 ! Specify as H2 nn_rxn_type(1) = 3 ! Specify as NH3 nn_rxn_type(2) = 3 ! Specify as NH3 com_ins(1,1) = xm(numpart_n2,1,nph) ! save x,y,z com position of N2 com_ins(2,1) = ym(numpart_n2,1,nph) com_ins(3,1) = zm(numpart_n2,1,nph) com_ins(1,2) = xm(numpart_h2a,2,nph) ! save x,y,z com position of H2a com_ins(2,2) = ym(numpart_h2a,2,nph) com_ins(3,2) = zm(numpart_h2a,2,nph) ! generate an initial orientation for the inserted NH3 molecules ! q_ins(x/y/z,point_chg#,nh3#) DO j=1,2 ! do for both molecules q_ins(1,1,j) = 1.0 ! Assign initial charge orientation for NH3 q_ins(2,1,j) = 0.0 q_ins(3,1,j) = 0.0 q_ins(1,2,j) = COS(theta) q_ins(2,2,j) = SIN(theta)*COS(0.0) q_ins(3,2,j) = SIN(theta)*SIN(0.0) q_ins(1,3,j) = COS(theta) q_ins(2,3,j) = SIN(theta)*COS(0.667*pi) q_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) q_ins(1,4,j) = COS(theta) q_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) q_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) e_ins(1,1,j) = 1.0 ! Assign initial orientation for LJ sites on NH3 e_ins(2,1,j) = 0.0 e_ins(3,1,j) = 0.0 e_ins(1,2,j) = COS(theta) e_ins(2,2,j) = SIN(theta)*COS(0.0) e_ins(3,2,j) = SIN(theta)*SIN(0.0) e_ins(1,3,j) = COS(theta) e_ins(2,3,j) = SIN(theta)*COS(0.667*pi) e_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) e_ins(1,4,j) = COS(theta) e_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) e_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Give the NH3 molecules a random orientation. First set up the ! ! NH3 as in the initialization subroutine, then randomly rotate ! ! about each of the axes. ! ! ! ! orient1 will temporarily store the LJ orientations ! ! orient2 will temporarily store the charge orientations ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ntype = 3 ! specify the type to be ammonia DO j=1,nins_rxn ! loop over the two ammonia molecules nflips = 10 DO i=1,nflips Iaxis(i) = INT(ran3(iseed)*3.0) + 1 ! pick a random axis END DO axes: DO l=1,nflips SELECT CASE(Iaxis(l)) CASE(1) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around x-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = e_ins(1,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) + sinDg*e_ins(3,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) - sinDg*e_ins(2,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = q_ins(1,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) + sinDg*q_ins(3,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) - sinDg*q_ins(2,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(2) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around y-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) - sinDg*e_ins(3,i,j) orient1(2,i,j) = e_ins(2,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) + sinDg*e_ins(1,i,j) END DO DO i=1,Nc(ntype) ! rotate charges sites orient2(1,i,j) = cosDg*q_ins(1,i,j) - sinDg*q_ins(3,i,j) orient2(2,i,j) = q_ins(2,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) + sinDg*q_ins(1,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(3) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around z-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) + sinDg*e_ins(2,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) - sinDg*e_ins(1,i,j) orient1(3,i,j) = e_ins(3,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = cosDg*q_ins(1,i,j) + sinDg*q_ins(2,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) - sinDg*q_ins(1,i,j) orient2(3,i,j) = q_ins(3,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO END SELECT END DO axes END DO ! Send the information to the Rxn_Force subroutine to calculate the change in energy deltaU = Rxn_Force(nins_rxn,com_ins,e_ins,q_ins,ndel_rxn,n_rxn_part,n_rxn_type,nn_rxn_type,nph) IF(novr.EQ.1) RETURN ! return if there is an overlap in the configuration IF(deltaU/Temp(kk).GT.100) RETURN !! Acceptance criterion: The acceptance will be specific to the reaction IF (deltaU/Temp(kk).LE.-70.0) NU2accept = 1 IF (deltaU/Temp(kk).GT.-70.0) THEN Prs = REAL(Nb(1,nph))*REAL(Nb(2,nph))*REAL(Nb(2,nph)-1)*REAL(Nb(2,nph)-2) Prs = Prs/(REAL(Nb(3,nph)+2)*REAL(Nb(3,nph)+1)) Prs = Prs*EXP(-deltaU/Temp(kk)) Prs = Prs*partition_f(nph) z1 = MIN(1.0,Prs) z2 = ran3(iseed)*1.0 IF (z2.LT.z1) NU2accept = 1 END IF acceptance: IF(NU2accept.EQ.1) THEN ! Must re-order the lists to account for the forward reaction DO i=numpart_n2,Nb(1,nph)-1 xm(i,1,nph) = xm(i+1,1,nph) ! shift down the N2 x-com assignments ym(i,1,nph) = ym(i+1,1,nph) ! shift down the N2 y-com assignments zm(i,1,nph) = zm(i+1,1,nph) ! shift down the N2 z-com assignments DO j=1,Na(1) ! shift down the LJ site asignments for N2 ex(j,i,1,nph) = ex(j,i+1,1,nph) ey(j,i,1,nph) = ey(j,i+1,1,nph) ez(j,i,1,nph) = ez(j,i+1,1,nph) END DO DO j=1,Nc(1) ! shift down the charge site assignments qx(j,i,1,nph) = qx(j,i+1,1,nph) qy(j,i,1,nph) = qy(j,i+1,1,nph) qz(j,i,1,nph) = qz(j,i+1,1,nph) END DO END DO num_h2(1) = numpart_h2a ! temporarily store numbers for sorting num_h2(2) = numpart_h2b num_h2(3) = numpart_h2c IF(num_h2(2).GT.num_h2(3)) THEN ! Order the H2 assignment number num_h2(4) = num_h2(3) num_h2(3) = num_h2(2) num_h2(2) = num_h2(4) END IF IF(num_h2(1).GT.num_h2(2)) THEN num_h2(4) = num_h2(1) num_h2(1) = num_h2(2) num_h2(2) = num_h2(4) END IF IF(num_h2(2).GT.num_h2(3)) THEN num_h2(4) = num_h2(3) num_h2(3) = num_h2(2) num_h2(2) = num_h2(4) END IF number_1 = num_h2(1) number_2 = num_h2(2) number_3 = num_h2(3) DO i=number_1,number_2-2 xm(i,2,nph) = xm(i+1,2,nph) ! shift down the H2 x-com assignments ym(i,2,nph) = ym(i+1,2,nph) ! shift down the H2 y-com assignments zm(i,2,nph) = zm(i+1,2,nph) ! shift down the H2 z-com assignments DO j=1,Na(2) ! shift down the LJ site asignments for N2 ex(j,i,2,nph) = ex(j,i+1,2,nph) ey(j,i,2,nph) = ey(j,i+1,2,nph) ez(j,i,2,nph) = ez(j,i+1,2,nph) END DO DO j=1,Nc(2) ! shift down the charge site assignments qx(j,i,2,nph) = qx(j,i+1,2,nph) qy(j,i,2,nph) = qy(j,i+1,2,nph) qz(j,i,2,nph) = qz(j,i+1,2,nph) END DO END DO DO i=number_2-1,number_3-3 xm(i,2,nph) = xm(i+2,2,nph) ! shift down the H2 x-com assignments ym(i,2,nph) = ym(i+2,2,nph) ! shift down the H2 y-com assignments zm(i,2,nph) = zm(i+2,2,nph) ! shift down the H2 z-com assignments DO j=1,Na(2) ! shift down the LJ site asignments for N2 ex(j,i,2,nph) = ex(j,i+2,2,nph) ey(j,i,2,nph) = ey(j,i+2,2,nph) ez(j,i,2,nph) = ez(j,i+2,2,nph) END DO DO j=1,Nc(2) ! shift down the charge site assignments qx(j,i,2,nph) = qx(j,i+2,2,nph) qy(j,i,2,nph) = qy(j,i+2,2,nph) qz(j,i,2,nph) = qz(j,i+2,2,nph) END DO END DO DO i=number_3-2,Nb(2,nph)-3 xm(i,2,nph) = xm(i+3,2,nph) ! shift down the H2 x-com assignments ym(i,2,nph) = ym(i+3,2,nph) ! shift down the H2 y-com assignments zm(i,2,nph) = zm(i+3,2,nph) ! shift down the H2 z-com assignments DO j=1,Na(2) ! shift down the LJ site asignments for N2 ex(j,i,2,nph) = ex(j,i+3,2,nph) ey(j,i,2,nph) = ey(j,i+3,2,nph) ez(j,i,2,nph) = ez(j,i+3,2,nph) END DO DO j=1,Nc(2) ! shift down the charge site assignments qx(j,i,2,nph) = qx(j,i+3,2,nph) qy(j,i,2,nph) = qy(j,i+3,2,nph) qz(j,i,2,nph) = qz(j,i+3,2,nph) END DO END DO ! Must now assign the two new NH3's xm(Nb(3,nph)+1,3,nph) = com_ins(1,1) ! NH3 #1 ym(Nb(3,nph)+1,3,nph) = com_ins(2,1) zm(Nb(3,nph)+1,3,nph) = com_ins(3,1) xm(Nb(3,nph)+2,3,nph) = com_ins(1,2) ! NH3 #2 ym(Nb(3,nph)+2,3,nph) = com_ins(2,2) zm(Nb(3,nph)+2,3,nph) = com_ins(3,2) DO i=1,Na(3) ! assign the new LJ sites ex(i,Nb(3,nph)+1,3,nph) = e_ins(1,i,1) ! NH3 #1 ey(i,Nb(3,nph)+1,3,nph) = e_ins(2,i,1) ez(i,Nb(3,nph)+1,3,nph) = e_ins(3,i,1) ex(i,Nb(3,nph)+2,3,nph) = e_ins(1,i,2) ! NH3 #2 ey(i,Nb(3,nph)+2,3,nph) = e_ins(2,i,2) ez(i,Nb(3,nph)+2,3,nph) = e_ins(3,i,2) END DO DO i=1,Nc(3) ! assign the new charge sites qx(i,Nb(3,nph)+1,3,nph) = q_ins(1,i,1) ! NH3 #1 qy(i,Nb(3,nph)+1,3,nph) = q_ins(2,i,1) qz(i,Nb(3,nph)+1,3,nph) = q_ins(3,i,1) qx(i,Nb(3,nph)+2,3,nph) = q_ins(1,i,2) ! NH3 #2 qy(i,Nb(3,nph)+2,3,nph) = q_ins(2,i,2) qz(i,Nb(3,nph)+2,3,nph) = q_ins(3,i,2) END DO Nb(1,nph) = Nb(1,nph) - 1 ! One N2 molecule has been consumed Nb(2,nph) = Nb(2,nph) - 3 ! Three H2 molecules have been consumed Nb(3,nph) = Nb(3,nph) + 2 ! One NH3 molecule has been created Uold(nph) = Uold(nph) + deltaU ! update energy IF (t.GT.t_equil) acc1 = acc1 + 1.0 ! update acceptance count END IF acceptance END SUBROUTINE Forward !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Reverse !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This subroutine will perform a reverse reaction step. For the ammonia ! ! synthesis reaction, one NH3 will be changed to N2 and the other NH3 ! ! will be changed to an H2. In addition, two molecules of H2 will be ! ! inserted randomly into the fluid. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION com_ins(3,4),q_ins(3,ncgs,4),e_ins(3,natoms,4),n_rxn_part(2),n_rxn_type(2) DIMENSION orient1(3,natoms,2),orient2(3,ncgs,2),nn_rxn_type(4) NU3accept=0 IF (t.GT.t_equil) att2=att2+1.0 If(Nb(3,nph).LT.2) RETURN ! Return if not enough ammonia molecules are present numpart_nh3a = INT(AINT(ran3(iseed)*REAL(Nb(3,nph)))) + 1 ! pick an NH3 at random 1 numpart_nh3b = INT(AINT(ran3(iseed)*REAL(Nb(3,nph)))) + 1 ! pick an NH3 at random IF(numpart_nh3b.EQ.numpart_nh3a) GOTO 1 ! make sure the molecules aren't the same nins_rxn = 4; ndel_rxn = 2 ! insert: 1 N2 and 3H2, delete 2 NH3 n_rxn_part(1) = numpart_nh3a ! save the NH3 molecules #'s n_rxn_part(2) = numpart_nh3b ! save the NH3 molecules #'s n_rxn_type(1) = 3 ! specify as NH3 (reactant) n_rxn_type(2) = 3 ! specify as NH3 (reactant) nn_rxn_type(1) = 1 ! specify as N2 (product) nn_rxn_type(2) = 2 ! specify as H2 (product) nn_rxn_type(3) = 2 ! specify as H2 (product) nn_rxn_type(4) = 2 ! specify as H2 (product) com_ins(1,1) = xm(numpart_nh3a,3,nph) com_ins(2,1) = ym(numpart_nh3a,3,nph) com_ins(3,1) = zm(numpart_nh3a,3,nph) com_ins(1,2) = xm(numpart_nh3b,3,nph) com_ins(2,2) = ym(numpart_nh3b,3,nph) com_ins(3,2) = zm(numpart_nh3b,3,nph) com_ins(1,3) = (ran3(iseed)-0.5)*box(1,nph) ! random pos for ins #3 com_ins(2,3) = (ran3(iseed)-0.5)*box(2,nph) com_ins(3,3) = (ran3(iseed)-0.5)*box(3,nph) com_ins(1,4) = (ran3(iseed)-0.5)*box(1,nph) ! random pos for ins #4 com_ins(2,4) = (ran3(iseed)-0.5)*box(2,nph) com_ins(3,4) = (ran3(iseed)-0.5)*box(3,nph) !! Generate initial orientations for the inserted N2 and H2 molecules ! DO i=1,4 q_ins(1,1,i) = 0.0 ! The 1st charge will be at the com q_ins(2,1,i) = 0.0 q_ins(3,1,i) = 0.0 ransq=2.0 ! This generates a random vector on a unit sphere DO WHILE (ransq.GE.1.0) ran1 = 1.0-2.0*ran3(iseed) ran2 = 1.0-2.0*ran3(iseed) ransq = ran1*ran1 + ran2*ran2 END DO ranh=2.0*SQRT(1.0-ransq) exnew=ran1*ranh eynew=ran2*ranh eznew=(1.0-2.0*ransq) q_ins(1,2,i) = exnew ! orientation for 2nd charge q_ins(2,2,i) = eynew q_ins(3,2,i) = eznew q_ins(1,3,i) = -exnew ! charges are symmetric on N2 and H2 q_ins(2,3,i) = -eynew q_ins(3,3,i) = -eznew e_ins(1,1,i) = exnew ! LJ site #1 e_ins(2,1,i) = eynew e_ins(3,1,i) = eznew e_ins(1,2,i) = -exnew ! LJ site #2 e_ins(2,2,i) = -eynew e_ins(3,2,i) = -eznew END DO deltaU = Rxn_Force(nins_rxn,com_ins,e_ins,q_ins,ndel_rxn,n_rxn_part,n_rxn_type,nn_rxn_type,nph) IF(novr.EQ.1) RETURN ! return if there is an overlap in the configuration IF(deltaU/Temp(kk).GT.100) RETURN !! Acceptance criterion: The acceptance will be specific to the reaction IF (deltaU/Temp(kk).LE.-70.0) NU3accept = 1 IF (deltaU/Temp(kk).GT.-70.0) THEN Prs = REAL(Nb(3,nph))*REAL(Nb(3,nph)-1) Prs = Prs/(REAL(Nb(1,nph)+1)*REAL(Nb(2,nph)+3)*REAL(Nb(2,nph)+2)*REAL(Nb(2,nph)+1)) Prs = Prs*EXP(-deltaU/Temp(kk)) Prs = Prs*partition_r(nph) z1 = MIN(1.0,Prs) z2 = ran3(iseed)*1.0 IF (z2.LT.z1) NU3accept = 1 END IF acceptance: IF(NU3accept.EQ.1) THEN number_1 = MIN(numpart_nh3a,numpart_nh3b) number_2 = MAX(numpart_nh3a,numpart_nh3b) nntyp = 3 ! this is for shifting ammonia assignments DO i=number_1,number_2-2 xm(i,nntyp,nph) = xm(i+1,nntyp,nph) ! shift down the NH3 x-com assignments ym(i,nntyp,nph) = ym(i+1,nntyp,nph) ! shift down the NH3 y-com assignments zm(i,nntyp,nph) = zm(i+1,nntyp,nph) ! shift down the NH3 z-com assignments DO j=1,Na(nntyp) ! shift down the LJ site asignments for NH3 ex(j,i,nntyp,nph) = ex(j,i+1,nntyp,nph) ey(j,i,nntyp,nph) = ey(j,i+1,nntyp,nph) ez(j,i,nntyp,nph) = ez(j,i+1,nntyp,nph) END DO DO j=1,Nc(nntyp) ! shift down the charge site assignments qx(j,i,nntyp,nph) = qx(j,i+1,nntyp,nph) qy(j,i,nntyp,nph) = qy(j,i+1,nntyp,nph) qz(j,i,nntyp,nph) = qz(j,i+1,nntyp,nph) END DO END DO DO i=number_2-1,Nb(nntyp,nph)-2 xm(i,nntyp,nph) = xm(i+2,nntyp,nph) ! shift down the NH3 x-com assignments ym(i,nntyp,nph) = ym(i+2,nntyp,nph) ! shift down the NH3 y-com assignments zm(i,nntyp,nph) = zm(i+2,nntyp,nph) ! shift down the NH3 z-com assignments DO j=1,Na(nntyp) ! shift down the LJ site asignments for NH3 ex(j,i,nntyp,nph) = ex(j,i+2,nntyp,nph) ey(j,i,nntyp,nph) = ey(j,i+2,nntyp,nph) ez(j,i,nntyp,nph) = ez(j,i+2,nntyp,nph) END DO DO j=1,Nc(nntyp) ! shift down the charge site assignments qx(j,i,nntyp,nph) = qx(j,i+2,nntyp,nph) qy(j,i,nntyp,nph) = qy(j,i+2,nntyp,nph) qz(j,i,nntyp,nph) = qz(j,i+2,nntyp,nph) END DO END DO ! Must now assign the new N2 xm(Nb(1,nph)+1,1,nph) = com_ins(1,1) ! new N2 com ym(Nb(1,nph)+1,1,nph) = com_ins(2,1) zm(Nb(1,nph)+1,1,nph) = com_ins(3,1) xm(Nb(2,nph)+1,2,nph) = com_ins(1,2) ! new H2 #1 com ym(Nb(2,nph)+1,2,nph) = com_ins(2,2) zm(Nb(2,nph)+1,2,nph) = com_ins(3,2) xm(Nb(2,nph)+2,2,nph) = com_ins(1,3) ! new H2 #2 com ym(Nb(2,nph)+2,2,nph) = com_ins(2,3) zm(Nb(2,nph)+2,2,nph) = com_ins(3,3) xm(Nb(2,nph)+3,2,nph) = com_ins(1,4) ! new H2 #3 com ym(Nb(2,nph)+3,2,nph) = com_ins(2,4) zm(Nb(2,nph)+3,2,nph) = com_ins(3,4) DO i=1,Na(1) ! assign the new LJ sites ex(i,Nb(1,nph)+1,1,nph) = e_ins(1,i,1) ! N2 ey(i,Nb(1,nph)+1,1,nph) = e_ins(2,i,1) ez(i,Nb(1,nph)+1,1,nph) = e_ins(3,i,1) END DO DO i=1,Na(2) ! assign the new LJ sites ex(i,Nb(2,nph)+1,2,nph) = e_ins(1,i,2) ! H2 #1 ey(i,Nb(2,nph)+1,2,nph) = e_ins(2,i,2) ez(i,Nb(2,nph)+1,2,nph) = e_ins(3,i,2) ex(i,Nb(2,nph)+2,2,nph) = e_ins(1,i,3) ! H2 #2 ey(i,Nb(2,nph)+2,2,nph) = e_ins(2,i,3) ez(i,Nb(2,nph)+2,2,nph) = e_ins(3,i,3) ex(i,Nb(2,nph)+3,2,nph) = e_ins(1,i,4) ! H2 #3 ey(i,Nb(2,nph)+3,2,nph) = e_ins(2,i,4) ez(i,Nb(2,nph)+3,2,nph) = e_ins(3,i,4) END DO DO i=1,Nc(1) ! assign the new charge sites qx(i,Nb(1,nph)+1,1,nph) = q_ins(1,i,1) ! N2 qy(i,Nb(1,nph)+1,1,nph) = q_ins(2,i,1) qz(i,Nb(1,nph)+1,1,nph) = q_ins(3,i,1) END DO DO i=1,Nc(2) ! assign the new LJ sites qx(i,Nb(2,nph)+1,2,nph) = q_ins(1,i,2) ! H2 #1 qy(i,Nb(2,nph)+1,2,nph) = q_ins(2,i,2) qz(i,Nb(2,nph)+1,2,nph) = q_ins(3,i,2) qx(i,Nb(2,nph)+2,2,nph) = q_ins(1,i,3) ! H2 #2 qy(i,Nb(2,nph)+2,2,nph) = q_ins(2,i,3) qz(i,Nb(2,nph)+2,2,nph) = q_ins(3,i,3) qx(i,Nb(2,nph)+3,2,nph) = q_ins(1,i,4) ! H2 #3 qy(i,Nb(2,nph)+3,2,nph) = q_ins(2,i,4) qz(i,Nb(2,nph)+3,2,nph) = q_ins(3,i,4) END DO Nb(1,nph) = Nb(1,nph) + 1 ! One N2 molecule has been consumed Nb(2,nph) = Nb(2,nph) + 3 ! Three H2 molecules have been consumed Nb(3,nph) = Nb(3,nph) - 2 ! One NH3 molecule has been created Uold(nph) = Uold(nph) + deltaU ! update energy IF (t.GT.t_equil) acc2 = acc2 + 1.0 ! update acceptance counter END IF acceptance END SUBROUTINE Reverse !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Exchange_1_2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This subroutine will move a molecule from the bulk phase(1) to ! ! the pore phase(2). For the ammonia synthesis, only ammonia will need ! ! to be transfered, OR nitrogen and hydrogen may be transfered. I will ! ! start by swapping only ammonia. Later I will try swapping only ! ! nitrogen and hydrogen to see which method equilibrates faster. I will ! ! try to use the same force subroutine as the reaction moves. ! ! ! ! I believe (as a very late afterthought) that H2 and N2 need to be ! ! swapped, not just NH3. Otherwise the N2:H2 ratio will remain ! ! identical in both phases and will be equal to the initial ratio ! ! 1:3 as I have set up. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION com_ins(3,4),q_ins(3,ncgs,4),e_ins(3,natoms,4),n_rxn_part(4),n_rxn_type(4) DIMENSION orient1(3,natoms,2),orient2(3,ncgs,2),nn_rxn_type(4) DIMENSION com_ins2(3,4),q_ins2(3,ncgs,4),e_ins2(3,natoms,4),n_rxn_part2(4),n_rxn_type2(4) DIMENSION nn_rxn_type2(4),Iaxis(10) NU4accept=0;xold=0.0;yold=0.0;zold=0.0;nonvr=0.0;deltaU=0.0 xnew=0.0;ynew=0.0;znew=0.0;Vz=0.0;zrand=0.0;U_in=0.0;U_out=0.0 Vol(2)=box(1,2)*box(2,2)*box(3,2) Vol(1)=box(1,1)*box(2,1)*box(3,1) !!!!!!!!!!!!!!!!!!!!! RANDOM MOLECULE DELETION FROM BULK PHASE !!!!!!!!!!!!!!!!! ! n_rxn_type2(1) = INT(ran3(iseed)*3.0) + 1 ! randomly pick N2, H2, or NH3 to swap ! n_rxn_type2(1) = INT(ran3(iseed)*2.0) + 2 ! randomly pick H2 or NH3 to swap n_rxn_type2(1) = INT(ran3(iseed)*2.0) + 1 ! randomly pick N2 or H2 to swap ! n_rxn_type2(1) = 1 ! only swap N2 IF((n_rxn_type2(1).EQ.1).AND.(Nb(1,1).LT.1)) RETURN IF((n_rxn_type2(1).EQ.2).AND.(Nb(2,1).LT.1)) RETURN IF((n_rxn_type2(1).EQ.3).AND.(Nb(3,1).LT.1)) RETURN n_rxn_part2(1) = INT(AINT(ran3(iseed)*REAL(Nb(n_rxn_type2(1),1))))+1 ! pick exact molecule to delete nn_rxn_type2 = 0 ! insert 0 molecules ndel_rxn = 1 ! delete 1 molecule nins_rxn = 0 ! insert 0 molecules nph = 1 ! the particle will be deleted from the bulk com_ins2 = 0.0 ! zero all unnecessary parameters q_ins2 = 0.0 e_ins2 = 0.0 U_out = Rxn_Force(nins_rxn,com_ins2,e_ins2,q_ins2,ndel_rxn,n_rxn_part2,n_rxn_type2,nn_rxn_type2,nph) !!!!!!!!!!!!!!!!! RANDOM MOLECULE INSERTION INTO PORE PHASE !!!!!!!!!!!!!!!!!!! ! q_ins(x/y/z,point_chg#,mol#) com_ins(1,1) = (ran3(iseed)-0.5)*box(1,2) ! create a new position for insertion com_ins(2,1) = (ran3(iseed)-0.5)*box(2,2) com_ins(3,1) = (ran3(iseed)-0.5)*box(3,2) IF(n_rxn_type2(1).EQ.1) THEN ! If N2 molecule picked ransq=2.0 ! Generates a random vector on a unit sphere DO WHILE (ransq.GE.1.0) ran1 = 1.0-2.0*ran3(iseed) ran2 = 1.0-2.0*ran3(iseed) ransq = ran1*ran1 + ran2*ran2 END DO ranh=2.0*SQRT(1.0-ransq) exnew=ran1*ranh eynew=ran2*ranh eznew=(1.0-2.0*ransq) e_ins(1,1,1) = exnew ! orientation for inserted N2 molecule (atom #1) e_ins(2,1,1) = eynew e_ins(3,1,1) = eznew e_ins(1,2,1) = -exnew ! orientation for inserted N2 molecule (atom #2) e_ins(2,2,1) = -eynew e_ins(3,2,1) = -eznew q_ins(1,1,1) = 0.0 ! N2 charge located at center of mass q_ins(2,1,1) = 0.0 q_ins(3,1,1) = 0.0 q_ins(1,2,1) = exnew ! atom #1 of N2 molecule q_ins(2,2,1) = eynew q_ins(3,2,1) = eznew q_ins(1,3,1) = -exnew ! atom #2 of N2 molecule q_ins(2,3,1) = -eynew q_ins(3,3,1) = -eznew n_rxn_part(1) = 0 ! mole # deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 1 ! insert a N2 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 2 ! particle will be added to pore phase ELSE IF(n_rxn_type2(1).EQ.2) THEN ! If H2 molecule picked (1-LJ molecule) n_rxn_part(1) = 0 ! mol# deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 2 ! insert a H2 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 2 ! particle will be added to pore phase ELSE IF(n_rxn_type2(1).EQ.3) THEN ! If NH3 molecule picked (4-LJ molecule) j=1 ! insert (1) ammonia molecule q_ins(1,1,j) = 1.0 ! Assign initial charge orientation for NH3 q_ins(2,1,j) = 0.0 q_ins(3,1,j) = 0.0 q_ins(1,2,j) = COS(theta) q_ins(2,2,j) = SIN(theta)*COS(0.0) q_ins(3,2,j) = SIN(theta)*SIN(0.0) q_ins(1,3,j) = COS(theta) q_ins(2,3,j) = SIN(theta)*COS(0.667*pi) q_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) q_ins(1,4,j) = COS(theta) q_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) q_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) e_ins(1,1,j) = 1.0 ! Assign initial orientation for LJ sites on NH3 e_ins(2,1,j) = 0.0 e_ins(3,1,j) = 0.0 e_ins(1,2,j) = COS(theta) e_ins(2,2,j) = SIN(theta)*COS(0.0) e_ins(3,2,j) = SIN(theta)*SIN(0.0) e_ins(1,3,j) = COS(theta) e_ins(2,3,j) = SIN(theta)*COS(0.667*pi) e_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) e_ins(1,4,j) = COS(theta) e_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) e_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Give the NH3 molecules a random orientation. First set up the ! ! NH3 as in the initialization subroutine, then randomly rotate ! ! about each of the axes. ! ! ! ! orient1 will temporarily store the LJ orientations ! ! orient2 will temporarily store the charge orientations ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ntype = 3 ! specify the type to be ammonia nflips = 10 ! perform 10 random rotations before inserting DO i=1,nflips Iaxis(i) = INT(ran3(iseed)*3.0) + 1 ! pick a random axis END DO axes: DO l=1,nflips SELECT CASE(Iaxis(l)) CASE(1) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around x-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = e_ins(1,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) + sinDg*e_ins(3,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) - sinDg*e_ins(2,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = q_ins(1,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) + sinDg*q_ins(3,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) - sinDg*q_ins(2,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(2) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around y-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) - sinDg*e_ins(3,i,j) orient1(2,i,j) = e_ins(2,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) + sinDg*e_ins(1,i,j) END DO DO i=1,Nc(ntype) ! rotate charges sites orient2(1,i,j) = cosDg*q_ins(1,i,j) - sinDg*q_ins(3,i,j) orient2(2,i,j) = q_ins(2,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) + sinDg*q_ins(1,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(3) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around z-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) + sinDg*e_ins(2,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) - sinDg*e_ins(1,i,j) orient1(3,i,j) = e_ins(3,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = cosDg*q_ins(1,i,j) + sinDg*q_ins(2,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) - sinDg*q_ins(1,i,j) orient2(3,i,j) = q_ins(3,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO END SELECT END DO axes !!!!!!!!!!!!!!!!!!!!!! END AMMONIA ROTATIONS !!!!!!!!!!!!!!!!!!!!!!!! n_rxn_part(1) = 0 ! mol# deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 3 ! insert a NH3 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 2 ! particle will be added to pore phase END IF IF(t.GT.t_equil) cum_chem_pot(nn_rxn_type(1),2) = cum_chem_pot(nn_rxn_type(1),2) + 1.0 ! counter for mu U_in = Rxn_Force(nins_rxn,com_ins,e_ins,q_ins,ndel_rxn,n_rxn_part,n_rxn_type,nn_rxn_type,nph) IF(novr.EQ.1) RETURN ! return if there is an overlap in the configuration deltaU = U_in + U_out ! Change in energy IF (deltaU/Temp(kk).GT.50.0) RETURN ! Acceptance criterion IF(t.GT.t_equil) THEN chem_pot(nn_rxn_type(1),2) = chem_pot(nn_rxn_type(1),2) + (Vol(2)/REAL(Nb(nn_rxn_type(1),2)+1.0))*EXP(-U_in/Temp(kk)) END IF IF (deltaU/Temp(kk).LE.-70.0) NU4accept=1 IF (deltaU.GT.-70.0) THEN Prt = REAL(Nb(n_rxn_type2(1),1))*Vol(2)/(REAL(Nb(n_rxn_type2(1),2)+1)*Vol(1)) Prt = Prt*EXP(-(deltaU)/Temp(kk)) z1 = MIN(1.0,Prt) z2 = ran3(iseed)*1.0 IF (z2.LT.z1) NU4accept = 1 END IF !!!!!!!!!!!!!!!!!!!!!!! UPDATE MOLECULE ASSIGMENTS !!!!!!!!!!!!!!!!!!!!!!!! IF (NU4accept.EQ.1) THEN ! Update if accepted IF (t.GT.t_equil) nins = nins +1 number_1 =n_rxn_part2(1) nntyp = n_rxn_type2(1) ! molecule type assignments to shift nph = 1 ! Shift occurs in the bulk phase DO i=number_1,Nb(nntyp,nph)-1 ! Cycle through particles in bulk phase xm(i,nntyp,nph) = xm(i+1,nntyp,nph) ! shift down the x-com assignments ym(i,nntyp,nph) = ym(i+1,nntyp,nph) ! shift down the y-com assignments zm(i,nntyp,nph) = zm(i+1,nntyp,nph) ! shift down the z-com assignments DO j=1,Na(nntyp) ! shift down the LJ site asignments for NH3 ex(j,i,nntyp,nph) = ex(j,i+1,nntyp,nph) ey(j,i,nntyp,nph) = ey(j,i+1,nntyp,nph) ez(j,i,nntyp,nph) = ez(j,i+1,nntyp,nph) END DO DO j=1,Nc(nntyp) ! shift down the charge site assignments qx(j,i,nntyp,nph) = qx(j,i+1,nntyp,nph) qy(j,i,nntyp,nph) = qy(j,i+1,nntyp,nph) qz(j,i,nntyp,nph) = qz(j,i+1,nntyp,nph) END DO END DO ! Must now assign the newly inserted molecule nph = 2 ! This is for the pore phase xm(Nb(nntyp,nph)+1,nntyp,nph) = com_ins(1,1) ym(Nb(nntyp,nph)+1,nntyp,nph) = com_ins(2,1) zm(Nb(nntyp,nph)+1,nntyp,nph) = com_ins(3,1) DO i=1,Na(nntyp) ! assign the new LJ sites ex(i,Nb(nntyp,nph)+1,nntyp,nph) = e_ins(1,i,1) ey(i,Nb(nntyp,nph)+1,nntyp,nph) = e_ins(2,i,1) ez(i,Nb(nntyp,nph)+1,nntyp,nph) = e_ins(3,i,1) END DO DO i=1,Nc(nntyp) ! assign the new charge sites qx(i,Nb(nntyp,nph)+1,nntyp,nph) = q_ins(1,i,1) qy(i,Nb(nntyp,nph)+1,nntyp,nph) = q_ins(2,i,1) qz(i,Nb(nntyp,nph)+1,nntyp,nph) = q_ins(3,i,1) END DO Nb(nntyp,1) = Nb(nntyp,1) - 1 ! Remove the particle to Phase #1 Nb(nntyp,2) = Nb(nntyp,2) + 1 ! Add the particle to Phase #2 Uold(1) = Uold(1) + U_out ! particle energy removed from bulk Uold(2) = Uold(2) + U_in ! particle energy added to pore END IF END SUBROUTINE Exchange_1_2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Exchange_2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This subroutine will move a particle from the pore phase(2) to ! ! the bulk phase(1). For the ammonia synthesis, only ammonia will need ! ! to be transfered, OR nitrogen and hydrogen may be transfered. I will ! ! start by swapping only ammonia. Later I will try swapping only ! ! nitrogen and hydrogen to see which method equilibrates faster. I will ! ! try to use the same force subroutine as the reaction moves. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION com_ins(3,4),q_ins(3,ncgs,4),e_ins(3,natoms,4),n_rxn_part(4),n_rxn_type(4) DIMENSION orient1(3,natoms,2),orient2(3,ncgs,2),nn_rxn_type(4) DIMENSION com_ins2(3,4),q_ins2(3,ncgs,4),e_ins2(3,natoms,4),n_rxn_part2(4),n_rxn_type2(4) DIMENSION nn_rxn_type2(4),Iaxis(10) com_ins=0.0;q_ins=0.0;e_ins=0.0;n_rxn_part=0;n_rxn_type=0 orient1=0.0;orient2=0.0;nn_rxn_type=0;com_ins2=0.0;q_ins2=0.0 e_ins2=0.0;n_rxn_part2=0;n_rxn_type2=0;orient12=0.0;orient22=0.0;nn_rxn_type2=0 NU5accept=0;xold=0.0;yold=0.0;zold=0.0;DeltaU=0.0;U_in=0.0;U_out=0.0 xnew=0.0;ynew=0.0;znew=0.0;Vz=0.0;zrand=0.0;novr=0;Vx=0.0;Prt=0.0 Vol(2)=box(1,2)*box(2,2)*box(3,2) Vol(1)=box(1,1)*box(2,1)*box(3,1) !!!!!!!!!!!!!!!!!!!!! RANDOM MOLECULE DELETION FROM PORE PHASE !!!!!!!!!!!!!!!!! ! n_rxn_type2(1) = INT(ran3(iseed)*3.0) + 1 ! randomly pick N2, H2, or NH3 to swap n_rxn_type2(1) = INT(ran3(iseed)*2.0) + 1 ! randomly pick N2 or H2 to swap ! n_rxn_type2(1) = INT(ran3(iseed)*2.0) + 2 ! randomly pick H2 or NH3 ! n_rxn_type2(1) = 1 ! only swap N2 IF((n_rxn_type2(1).EQ.1).AND.(Nb(1,2).LT.1)) RETURN IF((n_rxn_type2(1).EQ.2).AND.(Nb(2,2).LT.1)) RETURN IF((n_rxn_type2(1).EQ.3).AND.(Nb(3,2).LT.1)) RETURN n_rxn_part2(1) = INT(AINT(ran3(iseed)*REAL(Nb(n_rxn_type2(1),2))))+1 ! pick exact molecule to delete nn_rxn_type2 = 0 ! insert 0 molecules ndel_rxn = 1 ! delete 1 molecule nins_rxn = 0 ! insert 0 molecules nph = 2 ! molecule will be deleted from pore phase com_ins2 = 0.0 ! zero all unnecessary parameters q_ins2 = 0.0 e_ins2 = 0.0 U_out = Rxn_Force(nins_rxn,com_ins2,e_ins2,q_ins2,ndel_rxn,n_rxn_part2,n_rxn_type2,nn_rxn_type2,nph) !!!!!!!!!!!!!!!!! RANDOM MOLECULE INSERTION INTO BULK PHASE !!!!!!!!!!!!!!!!!!! ! q_ins(x/y/z,point_chg#,nh3#) com_ins(1,1) = (ran3(iseed)-0.5)*box(1,1) ! create a new position in bulk phase com_ins(2,1) = (ran3(iseed)-0.5)*box(2,1) com_ins(3,1) = (ran3(iseed)-0.5)*box(3,1) IF(n_rxn_type2(1).EQ.1) THEN ! If N2 molecule picked ransq=2.0 ! Generates a random vector on a unit sphere DO WHILE (ransq.GE.1.0) ran1 = 1.0-2.0*ran3(iseed) ran2 = 1.0-2.0*ran3(iseed) ransq = ran1*ran1 + ran2*ran2 END DO ranh=2.0*SQRT(1.0-ransq) exnew=ran1*ranh eynew=ran2*ranh eznew=(1.0-2.0*ransq) e_ins(1,1,1) = exnew ! orientation for inserted N2 molecule (atom #1) e_ins(2,1,1) = eynew e_ins(3,1,1) = eznew e_ins(1,2,1) = -exnew ! orientation for inserted N2 molecule (atom #2) e_ins(2,2,1) = -eynew e_ins(3,2,1) = -eznew q_ins(1,1,1) = 0.0 ! N2 charge located at center of mass q_ins(2,1,1) = 0.0 q_ins(3,1,1) = 0.0 q_ins(1,2,1) = exnew ! atom #1 of N2 molecule q_ins(2,2,1) = eynew q_ins(3,2,1) = eznew q_ins(1,3,1) = -exnew ! atom #2 of N2 molecule q_ins(2,3,1) = -eynew q_ins(3,3,1) = -eznew n_rxn_part(1) = 0 ! mole # deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 1 ! insert a N2 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 1 ! particle will be added to bulk phase ELSE IF(n_rxn_type2(1).EQ.2) THEN ! If H2 molecule picked (1-LJ molecule) n_rxn_part(1) = 0 ! mol# deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 2 ! insert a H2 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 1 ! particle will be added to bulk phase ELSE IF(n_rxn_type2(1).EQ.3) THEN j=1 ! insert (1) ammonia molecule q_ins(1,1,j) = 1.0 ! Assign initial charge orientation for NH3 q_ins(2,1,j) = 0.0 q_ins(3,1,j) = 0.0 q_ins(1,2,j) = COS(theta) q_ins(2,2,j) = SIN(theta)*COS(0.0) q_ins(3,2,j) = SIN(theta)*SIN(0.0) q_ins(1,3,j) = COS(theta) q_ins(2,3,j) = SIN(theta)*COS(0.667*pi) q_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) q_ins(1,4,j) = COS(theta) q_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) q_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) e_ins(1,1,j) = 1.0 ! Assign initial orientation for LJ sites on NH3 e_ins(2,1,j) = 0.0 e_ins(3,1,j) = 0.0 e_ins(1,2,j) = COS(theta) e_ins(2,2,j) = SIN(theta)*COS(0.0) e_ins(3,2,j) = SIN(theta)*SIN(0.0) e_ins(1,3,j) = COS(theta) e_ins(2,3,j) = SIN(theta)*COS(0.667*pi) e_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) e_ins(1,4,j) = COS(theta) e_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) e_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Give the NH3 molecules a random orientation. First set up the ! ! NH3 as in the initialization subroutine, then randomly rotate ! ! about each of the axes. ! ! ! ! orient1 will temporarily store the LJ orientations ! ! orient2 will temporarily store the charge orientations ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ntype = 3 ! specify the type to be ammonia nflips = 10 ! perform 10 random rotations before inserting DO i=1,nflips Iaxis(i) = INT(ran3(iseed)*3.0) + 1 ! pick a random axis END DO axes: DO l=1,nflips SELECT CASE(Iaxis(l)) CASE(1) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around x-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = e_ins(1,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) + sinDg*e_ins(3,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) - sinDg*e_ins(2,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = q_ins(1,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) + sinDg*q_ins(3,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) - sinDg*q_ins(2,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(2) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around y-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) - sinDg*e_ins(3,i,j) orient1(2,i,j) = e_ins(2,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) + sinDg*e_ins(1,i,j) END DO DO i=1,Nc(ntype) ! rotate charges sites orient2(1,i,j) = cosDg*q_ins(1,i,j) - sinDg*q_ins(3,i,j) orient2(2,i,j) = q_ins(2,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) + sinDg*q_ins(1,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(3) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around z-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) + sinDg*e_ins(2,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) - sinDg*e_ins(1,i,j) orient1(3,i,j) = e_ins(3,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = cosDg*q_ins(1,i,j) + sinDg*q_ins(2,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) - sinDg*q_ins(1,i,j) orient2(3,i,j) = q_ins(3,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO END SELECT END DO axes !!!!!!!!!!!!!!!!!!!!!! END AMMONIA ROTATIONS !!!!!!!!!!!!!!!!!!!!!!!! n_rxn_part(1) = 0 ! mol# deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 3 ! insert a type(3) NH3 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 1 ! particle will be added to bulk phase END IF IF(t.GT.t_equil) cum_chem_pot(nn_rxn_type(1),1) = cum_chem_pot(nn_rxn_type(1),1) + 1.0 ! counter for mu U_in = Rxn_Force(nins_rxn,com_ins,e_ins,q_ins,ndel_rxn,n_rxn_part,n_rxn_type,nn_rxn_type,nph) IF(novr.EQ.1) RETURN ! return if there is an overlap in the configuration !!!!!!!!!!!!!!!!!!!!!!!!!! ACCEPTANCE CRITERION !!!!!!!!!!!!!!!!!!!!!!!!!! deltaU = U_in + U_out ! Must calculate the change in energy IF (deltaU/Temp(kk).GT.50.0) RETURN ! Acceptance criterion IF(t.GT.t_equil) THEN chem_pot(nn_rxn_type(1),1) = chem_pot(nn_rxn_type(1),1) + (Vol(1)/REAL(Nb(nn_rxn_type(1),1)+1.0))*EXP(-U_in/Temp(kk)) END IF IF (deltaU/Temp(kk).LE.-70.0) NU5accept=1 IF (deltaU/Temp(kk).GT.-70.0) THEN Prt = REAL(Nb(n_rxn_type2(1),2))*Vol(1)/(REAL(Nb(n_rxn_type2(1),1)+1)*Vol(2)) Prt = Prt*EXP(-deltaU/Temp(kk)) IF (ran3(iseed).LT.(MIN(1.0,Prt))) NU5accept = 1 END IF !!!!!!!!!!!!!!!!!!!!!!!!!!! UPDATE MOLECULE ASSIGNMENTS !!!!!!!!!!!!!!!!!!!!!!!!! IF (NU5accept.EQ.1) THEN ! Update if accepted IF (t.GT.t_equil) ndel = ndel +1 number_1 =n_rxn_part2(1) nntyp = n_rxn_type2(1) ! molecule type assignments to shift nph = 2 ! Shift occurs in the pore phase DO i=number_1,Nb(nntyp,nph)-1 ! Cycle through particles in pore phase xm(i,nntyp,nph) = xm(i+1,nntyp,nph) ! shift down the x-com assignments ym(i,nntyp,nph) = ym(i+1,nntyp,nph) ! shift down the y-com assignments zm(i,nntyp,nph) = zm(i+1,nntyp,nph) ! shift down the z-com assignments DO j=1,Na(nntyp) ! shift down the LJ site asignments ex(j,i,nntyp,nph) = ex(j,i+1,nntyp,nph) ey(j,i,nntyp,nph) = ey(j,i+1,nntyp,nph) ez(j,i,nntyp,nph) = ez(j,i+1,nntyp,nph) END DO DO j=1,Nc(nntyp) ! shift down the charge site assignments qx(j,i,nntyp,nph) = qx(j,i+1,nntyp,nph) qy(j,i,nntyp,nph) = qy(j,i+1,nntyp,nph) qz(j,i,nntyp,nph) = qz(j,i+1,nntyp,nph) END DO END DO ! Must now assign the newly inserted molecule nph = 1 ! This is for the bulk phase xm(Nb(nntyp,nph)+1,nntyp,nph) = com_ins(1,1) ym(Nb(nntyp,nph)+1,nntyp,nph) = com_ins(2,1) zm(Nb(nntyp,nph)+1,nntyp,nph) = com_ins(3,1) DO i=1,Na(nntyp) ! assign the new LJ sites ex(i,Nb(nntyp,nph)+1,nntyp,nph) = e_ins(1,i,1) ey(i,Nb(nntyp,nph)+1,nntyp,nph) = e_ins(2,i,1) ez(i,Nb(nntyp,nph)+1,nntyp,nph) = e_ins(3,i,1) END DO DO i=1,Nc(nntyp) ! assign the new charge sites qx(i,Nb(nntyp,nph)+1,nntyp,nph) = q_ins(1,i,1) qy(i,Nb(nntyp,nph)+1,nntyp,nph) = q_ins(2,i,1) qz(i,Nb(nntyp,nph)+1,nntyp,nph) = q_ins(3,i,1) END DO Nb(nntyp,2) = Nb(nntyp,2) - 1 ! Remove the particle from Phase #2 Nb(nntyp,1) = Nb(nntyp,1) + 1 ! Add the particle to Phase #1 Uold(1) = Uold(1) + U_in Uold(2) = Uold(2) + U_out END IF END SUBROUTINE Exchange_2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Chemical_potential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This subroutine will calculate the chemical potential when only ONE ! ! phase is being simulated. Otherwise, the particle swaps between the ! ! two phases will be adequate to calculate the chemical potential, and ! ! this subroutine will be unnecessary. I will assume that the single ! ! phase being studied is the bulk (1) phase. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION com_ins(3,4),q_ins(3,4,4),e_ins(3,4,4),n_rxn_part(nspecies),n_rxn_type(nspecies) DIMENSION orient1(3,4,2),orient2(3,4,2),nn_rxn_type(nspecies) DIMENSION com_ins2(3,4),q_ins2(3,4,4),e_ins2(3,4,4),n_rxn_part2(nspecies),n_rxn_type2(nspecies) DIMENSION nn_rxn_type2(nspecies),Iaxis(10) com_ins=0.0;q_ins=0.0;e_ins=0.0;n_rxn_part=0;n_rxn_type=0 orient1=0.0;orient2=0.0;nn_rxn_type=0;com_ins2=0.0;q_ins2=0.0 e_ins2=0.0;n_rxn_part2=0;n_rxn_type2=0;orient12=0.0;orient22=0.0;nn_rxn_type2=0 NU5accept=0;xold=0.0;yold=0.0;zold=0.0;DeltaU=0.0;U_in=0.0;U_out=0.0 xnew=0.0;ynew=0.0;znew=0.0;Vz=0.0;zrand=0.0;novr=0;Vx=0.0;Prt=0.0 Vol(2)=box(1,2)*box(2,2)*box(3,2) Vol(1)=box(1,1)*box(2,1)*box(3,1) IF(t.LT.t_equil) RETURN !!!!!!!!!!!!!!!!! RANDOM TEST-PARTICLE INSERTION INTO BULK PHASE !!!!!!!!!!!!!!!!!!! n_rxn_type2(1) = INT(ran3(iseed)*3.0) + 1 ! randomly pick N2, H2, or NH3 to insert com_ins(1,1) = (ran3(iseed)-0.5)*box(1,1) ! create a new position in bulk phase com_ins(2,1) = (ran3(iseed)-0.5)*box(2,1) com_ins(3,1) = (ran3(iseed)-0.5)*box(3,1) IF(n_rxn_type2(1).EQ.1) THEN ! If N2 molecule picked ransq=2.0 ! Generates a random vector on a unit sphere DO WHILE (ransq.GE.1.0) ran1 = 1.0-2.0*ran3(iseed) ran2 = 1.0-2.0*ran3(iseed) ransq = ran1*ran1 + ran2*ran2 END DO ranh=2.0*SQRT(1.0-ransq) exnew=ran1*ranh eynew=ran2*ranh eznew=(1.0-2.0*ransq) e_ins(1,1,1) = exnew ! orientation for inserted N2 molecule (atom #1) e_ins(2,1,1) = eynew e_ins(3,1,1) = eznew e_ins(1,2,1) = -exnew ! orientation for inserted N2 molecule (atom #2) e_ins(2,2,1) = -eynew e_ins(3,2,1) = -eznew q_ins(1,1,1) = 0.0 ! N2 charge located at center of mass q_ins(2,1,1) = 0.0 q_ins(3,1,1) = 0.0 q_ins(1,2,1) = exnew ! atom #1 of N2 molecule q_ins(2,2,1) = eynew q_ins(3,2,1) = eznew q_ins(1,3,1) = -exnew ! atom #2 of N2 molecule q_ins(2,3,1) = -eynew q_ins(3,3,1) = -eznew n_rxn_part(1) = 0 ! mole # deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 1 ! insert a N2 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 1 ! particle will be added to bulk phase ELSE IF(n_rxn_type2(1).EQ.2) THEN ! If H2 molecule picked (1-LJ molecule) n_rxn_part(1) = 0 ! mol# deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 2 ! insert a H2 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 1 ! particle will be added to bulk phase ELSE IF(n_rxn_type2(1).EQ.3) THEN j=1 ! insert (1) ammonia molecule q_ins(1,1,j) = 1.0 ! Assign initial charge orientation for NH3 q_ins(2,1,j) = 0.0 q_ins(3,1,j) = 0.0 q_ins(1,2,j) = COS(theta) q_ins(2,2,j) = SIN(theta)*COS(0.0) q_ins(3,2,j) = SIN(theta)*SIN(0.0) q_ins(1,3,j) = COS(theta) q_ins(2,3,j) = SIN(theta)*COS(0.667*pi) q_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) q_ins(1,4,j) = COS(theta) q_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) q_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) e_ins(1,1,j) = 1.0 ! Assign initial orientation for LJ sites on NH3 e_ins(2,1,j) = 0.0 e_ins(3,1,j) = 0.0 e_ins(1,2,j) = COS(theta) e_ins(2,2,j) = SIN(theta)*COS(0.0) e_ins(3,2,j) = SIN(theta)*SIN(0.0) e_ins(1,3,j) = COS(theta) e_ins(2,3,j) = SIN(theta)*COS(0.667*pi) e_ins(3,3,j) = SIN(theta)*SIN(0.667*pi) e_ins(1,4,j) = COS(theta) e_ins(2,4,j) = SIN(theta)*COS(-0.667*pi) e_ins(3,4,j) = SIN(theta)*SIN(-0.667*pi) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Give the NH3 molecules a random orientation. First set up the ! ! NH3 as in the initialization subroutine, then randomly rotate ! ! about each of the axes. ! ! ! ! orient1 will temporarily store the LJ orientations ! ! orient2 will temporarily store the charge orientations ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ntype = 3 ! specify the type to be ammonia nflips = 10 ! perform 10 random rotations before inserting DO i=2,nflips Iaxis(i) = INT(ran3(iseed)*2.0) + 2 ! pick a random axis END DO axes: DO l=1,nflips SELECT CASE(Iaxis(l)) CASE(1) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around x-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = e_ins(1,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) + sinDg*e_ins(3,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) - sinDg*e_ins(2,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = q_ins(1,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) + sinDg*q_ins(3,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) - sinDg*q_ins(2,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(2) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around y-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) - sinDg*e_ins(3,i,j) orient1(2,i,j) = e_ins(2,i,j) orient1(3,i,j) = cosDg*e_ins(3,i,j) + sinDg*e_ins(1,i,j) END DO DO i=1,Nc(ntype) ! rotate charges sites orient2(1,i,j) = cosDg*q_ins(1,i,j) - sinDg*q_ins(3,i,j) orient2(2,i,j) = q_ins(2,i,j) orient2(3,i,j) = cosDg*q_ins(3,i,j) + sinDg*q_ins(1,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO CASE(3) Dgamma= 2.0*pi*(ran3(iseed) - 0.5) ! random rotation from -180 => +180 deg around z-axis cosDg = COS(Dgamma) sinDg = SIN(Dgamma) DO i=1,Na(ntype) ! rotate LJ sites orient1(1,i,j) = cosDg*e_ins(1,i,j) + sinDg*e_ins(2,i,j) orient1(2,i,j) = cosDg*e_ins(2,i,j) - sinDg*e_ins(1,i,j) orient1(3,i,j) = e_ins(3,i,j) END DO DO i=1,Nc(ntype) ! rotate charge sites orient2(1,i,j) = cosDg*q_ins(1,i,j) + sinDg*q_ins(2,i,j) orient2(2,i,j) = cosDg*q_ins(2,i,j) - sinDg*q_ins(1,i,j) orient2(3,i,j) = q_ins(3,i,j) END DO DO i=1,Na(ntype) ! normalize the LJ site orientations IF (ABS(bl(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient1(1,i,j)**2 + orient1(2,i,j)**2 + orient1(3,i,j)**2) e_ins(1,i,j) = orient1(1,i,j)*rnorm e_ins(2,i,j) = orient1(2,i,j)*rnorm e_ins(3,i,j) = orient1(3,i,j)*rnorm END DO DO i=1,Nc(ntype) ! normalize the charge orientations IF (ABS(ql(i,ntype)).LE.0.02) CYCLE rnorm = 1.0/SQRT(orient2(1,i,j)**2 + orient2(2,i,j)**2 + orient2(3,i,j)**2) q_ins(1,i,j) = orient2(1,i,j)*rnorm q_ins(2,i,j) = orient2(2,i,j)*rnorm q_ins(3,i,j) = orient2(3,i,j)*rnorm END DO END SELECT END DO axes !!!!!!!!!!!!!!!!!!!!!! END AMMONIA ROTATIONS !!!!!!!!!!!!!!!!!!!!!!!! n_rxn_part(1) = 0 ! mol# deleted n_rxn_type(1) = 0 ! # mols to be deleted nn_rxn_type(1)= 3 ! insert a type(3) NH3 molecule ndel_rxn = 0 ! delete 0 molecule nins_rxn = 1 ! insert 1 molecules nph = 1 ! particle will be added to bulk phase END IF cum_chem_pot(nn_rxn_type(1),1) = cum_chem_pot(nn_rxn_type(1),1) + 1.0 ! counter for mu U_in = Rxn_Force(nins_rxn,com_ins,e_ins,q_ins,ndel_rxn,n_rxn_part,n_rxn_type,nn_rxn_type,nph) IF(novr.EQ.1) RETURN ! return if there is an overlap in the configuration !!!!!!!!!!!!!!!!!!!!!!!!!! ACCUMULATE AVERAGES !!!!!!!!!!!!!!!!!!!!!!!!!! IF (U_in/Temp(kk).GT.50.0) RETURN ! Avoid overflow chem_pot(nn_rxn_type(1),1) = chem_pot(nn_rxn_type(1),1) + (Vol(1)/REAL(Nb(nn_rxn_type(1),1)+1.0))*EXP(-U_in/Temp(kk)) END SUBROUTINE Chemical_potential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Partition !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This subroutine will need to be called every time the temperature ! ! changes or the volume changes. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' ! Partition function for N2 Vol(nph) = box(1,nph)*box(2,nph)*box(3,nph) q(1,nph) = 226.35*(Temp(kk)**1.5)*Vol(nph) ! q(trans) q(1,nph) = q(1,nph)*6.3194*Temp(kk) ! q(rot) q(1,nph) = q(1,nph)/(1-EXP(-92.6923/Temp(kk))) ! q(vib) ! Partition function for H2 q(2,nph) = 4.3695*(Temp(kk)**1.5)*Vol(nph) ! q(trans) q(2,nph) = q(2,nph)*0.21336*Temp(kk) ! q(rot) q(2,nph) = q(2,nph)/(1-EXP(-170.742/Temp(kk))) ! q(vib) ! Partition function for NH3 q(3,nph) = 106.152*(Temp(kk)**1.5)*Vol(nph) ! q(trans) q(3,nph) = q(3,nph)*3.1944*(Temp(kk)**1.5) ! q(rot) q(3,nph) = q(3,nph)/((1-EXP(-131.868/Temp(kk)))*(1-EXP(-37.363/Temp(kk)))*(1-EXP(-134.066/Temp(kk)))& &*(1-EXP(-134.066/Temp(kk)))*(1-EXP(-64.011/Temp(kk)))*(1-EXP(-64.011/Temp(kk)))) ! q(vib) ! Now I will combine the electronic contributions (this avoids overflows in exp calculations) ! H-H: 103.2 kcal/mol, N-N: 225.1 kcal/mol, 3N-H: 276.8 kcal/mol (McQuarrie) partition_f(nph) = EXP(D(kk)/Temp(kk))*q(3,nph)*q(3,nph)/(q(1,nph)*q(2,nph)*q(2,nph)*q(2,nph)) partition_r(nph) = EXP(-D(kk)/Temp(kk))*q(1,nph)*q(2,nph)*q(2,nph)*q(2,nph)/(q(3,nph)*q(3,nph)) END SUBROUTINE Partition !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Boundary INCLUDE 'ammonia.inc' phases: DO k=1,2 ! loop over phases b_c = REAL(2-k) species: DO i=1,nspecies ! loop over all species molecule: DO j=1,Nb(i,k) ! loop ocer all molecules of species(i) xm(j,i,k) = xm(j,i,k) - box(1,k)*REAL(ANINT(xm(j,i,k)/box(1,k))) ym(j,i,k) = ym(j,i,k) - box(2,k)*REAL(ANINT(ym(j,i,k)/box(2,k))) zm(j,i,k) = zm(j,i,k) - box(3,k)*REAL(ANINT(zm(j,i,k)/box(3,k)))*b_c END DO molecule END DO species END DO phases END SUBROUTINE Boundary !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION Force_calc(xmp,ymp,zmp,exn,eyn,ezn,qxn,qyn,qzn,np) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This is the main force subroutine which calculates the total force ! ! over all LJ sites and all partial charges. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION :: xmp(max2,nspecies,N),ymp(max2,nspecies,N),zmp(max2,nspecies,N) DIMENSION :: exn(natoms,max2,nspecies,N),eyn(natoms,max2,nspecies,N) DIMENSION :: ezn(natoms,max2,nspecies,N),qxn(ncgs,max2,nspecies,N) DIMENSION :: qyn(ncgs,max2,nspecies,N),qzn(ncgs,max2,nspecies,N) novr=0 Corr=0.0; VirLJ=0.0;VirE=0.0;force_calc=0.0; St=0.0; h=0.0 Vatt=0.0; Vrep=0.0; Vcorr=0.0; E_pot=0.0 sr2=0.0; sr6=0.0; sr12=0.0 xai=0.0; yai=0.0; zai=0.0 xij=0.0; yij=0.0; zij=0.0 rijsq=0.0;Vir_a=0.0;Vir_b=0.0;Vir_c=0.0 nTot_mol_nph=0;St=0.0 nbin=0;radcut=25.0 ! parameters for N-N RDF DO i=1,nspecies nTot_mol_nph = nTot_mol_nph + Nb(i,np) END DO IF (nTot_mol_nph.LE.0) RETURN ! check for particles boxx=box(1,np); boxy=box(2,np); boxz=box(3,np) b_c = REAL(2-np) rho_m(np) = REAL(nTot_mol_nph)/(boxx*boxy*boxz) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! type_1: DO i=1,nspecies ! outer loop over species type type_2: DO j=i,nspecies ! loop over species (i) and higher check: IF(j.EQ.i) THEN ! don't over count if same type mol_1a: DO k=1,Nb(i,np)-1 xmk = xmp(k,i,np) ! mol_1 COM ymk = ymp(k,i,np) zmk = zmp(k,i,np) mol_2a: DO l=k+1,Nb(j,np) xml = xmp(l,j,np) yml = ymp(l,j,np) zml = zmp(l,j,np) xr = xmk - xml ! mol_1 to mol_2 COM distance yr = ymk - yml zr = zmk - zml xr = xr - boxx*ANINT(xr/boxx) ! Minimum image yr = yr - boxy*ANINT(yr/boxy) zr = zr - boxz*ANINT(zr/boxz)*b_c ! min. image in pore z-direction removed r2 = (xr*xr) + (yr*yr) + (zr*zr) IF (r2.LT.radcut) THEN ! RDF calculation for N-N nbin = INT(SQRT(r2)/dr)+1 IF(nbin.LE.maxbin) THEN histA(nbin,1) = histA(nbin,1) + 2.0 END IF END IF IF(r2.LT.0.10) THEN ! return if an overlap occurs novr=1 PRINT *,'overlap in main force!' RETURN END IF IF (r2.LT.el_sqcut) THEN IF (r2.LT.sqcut) THEN ! calculate site forces site_1a: DO ns_1=1,Na(i) ! atomic sites on molecule type (i) site_2a: DO ns_2=1,Na(j) ! atomic sites on molecule type (j) xr_1 = xr + bl(ns_1,i)*exn(ns_1,k,i,np) - bl(ns_2,j)*exn(ns_2,l,j,np) ! site-site distances yr_1 = yr + bl(ns_1,i)*eyn(ns_1,k,i,np) - bl(ns_2,j)*eyn(ns_2,l,j,np) zr_1 = zr + bl(ns_1,i)*ezn(ns_1,k,i,np) - bl(ns_2,j)*ezn(ns_2,l,j,np) r2_1 = (xr_1*xr_1) + (yr_1*yr_1) + (zr_1*zr_1) r2i_1 = (1/r2_1) r6i_1 = r2i_1*r2i_1*r2i_1 sig_ij_6 = sig_m(ns_1,i,ns_2,j)**6 Vatt = Vatt - eps_m(ns_1,i,ns_2,j)*r6i_1*sig_ij_6 ! update attractive LJ potential Vrep = Vrep + eps_m(ns_1,i,ns_2,j)*r6i_1*r6i_1*(sig_ij_6*sig_ij_6) ! update repulsive LJ potential VirLJ = VirLJ + eps_m(ns_1,i,ns_2,j)*(xr_1*xr + yr_1*yr + zr_1*zr)*(r6i_1*sig_ij_6)*(r6i_1*(sig_ij_6)-0.5)*r2i_1 END DO site_2a END DO site_1a END IF charge_1a: DO nc_1=1,Nc(i) ! charge sites on molecule type (i) charge_2a: DO nc_2=1,Nc(j) ! charge sites on molecule type (j) xr_1 = xr + ql(nc_1,i)*qxn(nc_1,k,i,np) - ql(nc_2,j)*qxn(nc_2,l,j,np) ! site-site distances yr_1 = yr + ql(nc_1,i)*qyn(nc_1,k,i,np) - ql(nc_2,j)*qyn(nc_2,l,j,np) zr_1 = zr + ql(nc_1,i)*qzn(nc_1,k,i,np) - ql(nc_2,j)*qzn(nc_2,l,j,np) r_1 = SQRT((xr_1*xr_1) + (yr_1*yr_1) + (zr_1*zr_1)) E_pot = E_pot + charge(nc_1,i)*charge(nc_2,j)/r_1 VirE = VirE + charge(nc_1,i)*charge(nc_2,j)*(xr_1*xr + yr_1*yr + zr_1*zr)/(r_1*r_1*r_1) END DO charge_2a END DO charge_1a END IF END DO mol_2a IF (np.EQ.2) THEN ! Calculate Molecule-Pore potential site_1a_st1: DO ns_1=1,Na(i) ! Only account for LJ sites zai = zmk + bl(ns_1,i)*ezn(ns_1,k,i,np) ! Position of each LJ site h = 0.5*box(3,2) ! Wall positions St = St + St_eps(ns_1,i)*St_sig(ns_1,i)*St_sig(ns_1,i)*(0.2*(St_sig(ns_1,i)/(h - zai))**10& & + 0.2*(St_sig(ns_1,i)/(h + zai))**10& & - 0.5*(St_sig(ns_1,i)/(h - zai))**4& & - 0.5*(St_sig(ns_1,i)/(h + zai))**4& & - St_sig(ns_1,i)**4/(6.0*delta*(h - zai + (0.61*delta))**3)& & - St_sig(ns_1,i)**4/(6.0*delta*(h + zai + (0.61*delta))**3)) END DO site_1a_st1 END IF END DO mol_1a IF (np.EQ.2) THEN ! Account for last Molecule-Pore potential IF(Nb(i,np).GE.1) THEN k=Nb(i,np) ! last molecule of type (i) site_1a_st2: DO ns_1=1,Na(i) ! Only account for LJ sites zai = zmp(k,i,np) + bl(ns_1,i)*ezn(ns_1,k,i,np) h = 0.5*box(3,2) St = St + St_eps(ns_1,i)*St_sig(ns_1,i)*St_sig(ns_1,i)*(0.2*(St_sig(ns_1,i)/(h - zai))**10& & + 0.2*(St_sig(ns_1,i)/(h + zai))**10& & - 0.5*(St_sig(ns_1,i)/(h - zai))**4& & - 0.5*(St_sig(ns_1,i)/(h + zai))**4& & - St_sig(ns_1,i)**4/(6.0*delta*(h - zai + (0.61*delta))**3)& & - St_sig(ns_1,i)**4/(6.0*delta*(h + zai + (0.61*delta))**3)) END DO site_1a_st2 END IF END IF ELSE ! loop over all mol if diff type mol_1b: DO k=1,Nb(i,np) xmk = xmp(k,i,np) ! mol_1 COM ymk = ymp(k,i,np) zmk = zmp(k,i,np) mol_2b: DO l=1,Nb(j,np) xml = xmp(l,j,np) ! mol_2 COM yml = ymp(l,j,np) zml = zmp(l,j,np) xr = xmk - xml ! mol_1 to mol_2 COM distance yr = ymk - yml zr = zmk - zml xr = xr - boxx*ANINT(xr/boxx) ! Minimum image yr = yr - boxy*ANINT(yr/boxy) zr = zr - boxz*ANINT(zr/boxz)*b_c ! min. image in pore z-direction removed r2 = (xr*xr) + (yr*yr) + (zr*zr) IF (r2.LT.radcut) THEN ! RDF calculation for N-N nbin = INT(SQRT(r2)/dr)+1 IF(nbin.LE.maxbin) THEN histA(nbin,1) = histA(nbin,1) + 2.0 END IF END IF IF(r2.LT.0.10) THEN ! return if an overlap occurs novr=1 PRINT *,'overlap in main force!' RETURN END IF IF (r2.LT.el_sqcut) THEN IF (r2.LT.sqcut) THEN ! calculate site forces site_1b: DO ns_1=1,Na(i) ! atomic sites on molecule type (i) site_2b: DO ns_2=1,Na(j) ! atomic sites on molecule type (j) xr_1 = xr + bl(ns_1,i)*exn(ns_1,k,i,np) - bl(ns_2,j)*exn(ns_2,l,j,np) ! site-site distances yr_1 = yr + bl(ns_1,i)*eyn(ns_1,k,i,np) - bl(ns_2,j)*eyn(ns_2,l,j,np) zr_1 = zr + bl(ns_1,i)*ezn(ns_1,k,i,np) - bl(ns_2,j)*ezn(ns_2,l,j,np) r2_1 = (xr_1*xr_1) + (yr_1*yr_1) + (zr_1*zr_1) r2i_1 = (1/r2_1) r6i_1 = r2i_1*r2i_1*r2i_1 sig_ij_6 = sig_m(ns_1,i,ns_2,j)**6 Vatt = Vatt - eps_m(ns_1,i,ns_2,j)*r6i_1*sig_ij_6 ! update attractive LJ potential Vrep = Vrep + eps_m(ns_1,i,ns_2,j)*r6i_1*r6i_1*(sig_ij_6*sig_ij_6) ! update repulsive LJ potential VirLJ = VirLJ + eps_m(ns_1,i,ns_2,j)*(xr_1*xr + yr_1*yr + zr_1*zr)*(r6i_1*sig_ij_6)*((r6i_1*sig_ij_6)-0.5)*r2i_1 END DO site_2b END DO site_1b END IF charge_1b: DO nc_1=1,Nc(i) ! charge sites on molecule type (i) charge_2b: DO nc_2=1,Nc(j) ! charge sites on molecule type (j) xr_1 = xr + ql(nc_1,i)*qxn(nc_1,k,i,np) - ql(nc_2,j)*qxn(nc_2,l,j,np) ! site-site distances yr_1 = yr + ql(nc_1,i)*qyn(nc_1,k,i,np) - ql(nc_2,j)*qyn(nc_2,l,j,np) zr_1 = zr + ql(nc_1,i)*qzn(nc_1,k,i,np) - ql(nc_2,j)*qzn(nc_2,l,j,np) r_1 = SQRT((xr_1*xr_1) + (yr_1*yr_1) + (zr_1*zr_1)) E_pot = E_pot + charge(nc_1,i)*charge(nc_2,j)/(r_1) VirE = VirE + charge(nc_1,i)*charge(nc_2,j)*(xr_1*xr + yr_1*yr + zr_1*zr)/(r_1*r_1*r_1) END DO charge_2b END DO charge_1b END IF END DO mol_2b END DO mol_1b END IF check END DO type_2 END DO type_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! force_calc = 4.0*(Vrep + Vatt) + 1379.09*E_pot + St*4.0*pi*delta*4.172 VirNew = 48.0*VirLJ + 1379.09*VirE P_temp = rho_m(1)*Temp(kk) + VirNew/(3.0*box(1,1)*box(2,1)*box(3,1)) END FUNCTION Force_calc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FUNCTION Rxn_Force(ninsert,com_i,e_i,q_i,ndelete,ids,ityp,iityp,np) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! This function will return the value of the change in the potential ! ! energy upon a reaction step. It is important that all of the ! ! energy is counted - this includes interactions BETWEEN the newly ! ! inserted molecules, and BETWEEN the deleted molecules. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! ninsert: # of molecules to be inserted during reaction ! ! com_i: x,y,z array of the com of the inserted molecules ! ! e_i: x,y,z array of orientation of LJ sites on inserted molecules ! ! q_i: x,y,z array of orientation of charge sites on inserted mol ! ! ndelete: # of molecules deleted during the reaction ! ! ids: array of molecule id #'s from Nb(i,nph) which will be deleted ! ! ityp: array of molecules types which will be deleted ! ! iityp: array of molecule types which will be inserted ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INCLUDE 'ammonia.inc' DIMENSION com_i(3,4),e_i(3,natoms,4),q_i(3,ncgs,4),ids(4),ityp(4),iityp(4) novr=0;dV=0.0;Vanew=0.0;Vrnew=0.0;Vaold=0.0;Vrold=0.0 Vatt=0.0;Vrep=0.0;BForce=0.0;zst=0.0;St_ins=0.0;St_del=0.0 h=box(3,2)*0.5;Vir_old=0.0;Vir_new=0.0 E_pot_old=0.0;E_pot_new=0.0 b_c = REAL(2-np) boxx=box(1,np); boxy=box(2,np); boxz=box(3,np) ! Calculate the energy of the deleted particles (N2 + 3H2) with the other particles ! do ii=1,ndelete ! loop over the deleted particles (ii) xmk = xm(ids(ii),ityp(ii),np) ! COM of old particle ymk = ym(ids(ii),ityp(ii),np) zmk = zm(ids(ii),ityp(ii),np) type_2a: DO j=1,nspecies ! loop over species types mol_2a: DO l=1,Nb(j,np) ! loop over molecules xml = xm(l,j,np) yml = ym(l,j,np) zml = zm(l,j,np) DO jj=1,ii-1 ! skip loop if already counted previously IF((l.EQ.ids(jj)).AND.(j.EQ.ityp(jj))) GOTO 100 END DO IF((l.EQ.ids(ii)).AND.(j.EQ.ityp(ii))) GOTO 100 ! skip over same-particle calculation xr = xmk - xml ! mol_1 to mol_2 COM distance yr = ymk - yml zr = zmk - zml xr = xr - boxx*ANINT(xr/boxx) ! Minimum image yr = yr - boxy*ANINT(yr/boxy) zr = zr - boxz*ANINT(zr/boxz)*b_c ! min. image in pore z-direction removed r2 = (xr*xr) + (yr*yr) + (zr*zr) IF (r2.LT.el_sqcut) THEN IF (r2.LT.sqcut) THEN ! calculate site forces IF(r2.LT.0.20) THEN ! stop if an overlap is already present PRINT *,'Overlap Present!!' PRINT *,'xyz1',xmk,ymk,zmk PRINT *,'xyz2',xml,yml,zml PRINT *,'Phase',nph STOP END IF site_1a: DO ns_1=1,Na(ityp(ii)) ! LJ sites on molecule (i) site_2a: DO ns_2=1,Na(j) ! LJ sites on molecule (j) xr_1 = xr + bl(ns_1,ityp(ii))*ex(ns_1,ids(ii),ityp(ii),np) - bl(ns_2,j)*ex(ns_2,l,j,np) ! site-site distances yr_1 = yr + bl(ns_1,ityp(ii))*ey(ns_1,ids(ii),ityp(ii),np) - bl(ns_2,j)*ey(ns_2,l,j,np) zr_1 = zr + bl(ns_1,ityp(ii))*ez(ns_1,ids(ii),ityp(ii),np) - bl(ns_2,j)*ez(ns_2,l,j,np) r2_1 = (xr_1*xr_1) + (yr_1*yr_1) + (zr_1*zr_1) r2i_1 = (1/r2_1) r6i_1 = r2i_1*r2i_1*r2i_1 sig_ij_6 = sig_m(ns_1,ityp(ii),ns_2,j)**6 Vaold = Vaold - (eps_m(ns_1,ityp(ii),ns_2,j)*r6i_1*sig_ij_6) ! update attractive LJ pot Vrold = Vrold + (eps_m(ns_1,ityp(ii),ns_2,j)*r6i_1*r6i_1*(sig_ij_6*sig_ij_6)) ! update repulsive LJ pot END DO site_2a END DO site_1a END IF charge_1a: DO nc_1=1,Nc(ityp(ii)) ! charge sites on molecule type (i) charge_2a: DO nc_2=1,Nc(j) ! charge sites on molecule type (j) xr_1 = xr + ql(nc_1,ityp(ii))*qx(nc_1,ids(ii),ityp(ii),np) - ql(nc_2,j)*qx(nc_2,l,j,np) ! site-site distances yr_1 = yr + ql(nc_1,ityp(ii))*qy(nc_1,ids(ii),ityp(ii),np) - ql(nc_2,j)*qy(nc_2,l,j,np) zr_1 = zr + ql(nc_1,ityp(ii))*qz(nc_1,ids(ii),ityp(ii),np) - ql(nc_2,j)*qz(nc_2,l,j,np) r_1 = SQRT((xr_1*xr_1) + (yr_1*yr_1) + (zr_1*zr