From fe94cd4c40414972a7569a2b17e3cd5b3cad116d Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Mon, 29 Apr 2024 22:59:28 +0200 Subject: [PATCH] Fixed some bugs in the new PRNG (basically direct calls to prng.rng fail) Also created examples of atomswap and alchemical swaps, so that at least they get tested for running properly --- drivers/f90/LJ.f90 | 82 +++++++++++++ drivers/f90/driver.f90 | 15 ++- .../features/monte_carlo/alchemical/init.xyz | 8 ++ .../features/monte_carlo/alchemical/input.xml | 42 +++++++ .../features/monte_carlo/alchemical/run.sh | 16 +++ .../features/monte_carlo/atomswap/input.xml | 42 +++++++ .../features/monte_carlo/atomswap/ljmix.xyz | 1 + examples/features/monte_carlo/atomswap/run.sh | 16 +++ examples/init_files/ljmix.xyz | 110 ++++++++++++++++++ ipi/engine/motion/al6xxx_kmc.py | 2 +- ipi/engine/motion/alchemy.py | 8 +- ipi/engine/motion/atomswap.py | 8 +- ipi/utils/prng.py | 40 ++++++- 13 files changed, 373 insertions(+), 17 deletions(-) create mode 100644 examples/features/monte_carlo/alchemical/init.xyz create mode 100644 examples/features/monte_carlo/alchemical/input.xml create mode 100644 examples/features/monte_carlo/alchemical/run.sh create mode 100644 examples/features/monte_carlo/atomswap/input.xml create mode 120000 examples/features/monte_carlo/atomswap/ljmix.xyz create mode 100644 examples/features/monte_carlo/atomswap/run.sh create mode 100644 examples/init_files/ljmix.xyz diff --git a/drivers/f90/LJ.f90 b/drivers/f90/LJ.f90 index c60fbacc4..67343ef08 100644 --- a/drivers/f90/LJ.f90 +++ b/drivers/f90/LJ.f90 @@ -209,4 +209,86 @@ SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list, END SUBROUTINE + SUBROUTINE LJMix_getall(n_type2, rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial) + ! Calculates the LJ potential energy and virial and the forces + ! acting on all the atoms, using 0.4 sigma and 0.4 epsilon for the first n_type2 atoms + ! + ! Args: + ! n_type2: First atom are considered as type 2 + ! rc: The cut-off radius. + ! sigma: The LJ distance parameter. NOTE: x0 of the spring is set to be equal to sigma + ! eps: The LJ energy parameter. NOTE: stiffness of the spring is set to be equal: 36*2^(2/3) eps + ! natoms: The number of atoms in the system. + ! atoms: A vector holding all the atom positions. + ! cell_h: The simulation box cell vector matrix. + ! cell_ih: The inverse of the simulation box cell vector matrix. + ! index_list: A array giving the last index of n_list that + ! gives the neighbours of a given atom. + ! n_list: An array giving the indices of the atoms that neighbour + ! the atom determined by index_list. + ! pot: The total potential energy of the system. + ! forces: An array giving the forces acting on all the atoms. + ! virial: The virial tensor, not divided by the volume. + + DOUBLE PRECISION, INTENT(IN) :: rc + DOUBLE PRECISION, INTENT(IN) :: sigma + DOUBLE PRECISION, INTENT(IN) :: eps + INTEGER, INTENT(IN) :: natoms + DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(IN) :: atoms + DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h + DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih + INTEGER, DIMENSION(natoms), INTENT(IN) :: index_list + INTEGER, DIMENSION(natoms*(natoms-1)/2), INTENT(IN) :: n_list + DOUBLE PRECISION, INTENT(OUT) :: pot + DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(OUT) :: forces + DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: virial + + INTEGER i, j, k, l, start, n_type2 + DOUBLE PRECISION, DIMENSION(3) :: fij, rij + DOUBLE PRECISION sigmaij, epsij, r2, pot_ij, pot_lr, vir_lr, volume + + start = 1 + DO i = 1, natoms - 1 + ! Only loops over the neigbour list, not all the atoms. + DO j = start, index_list(i) + CALL vector_separation(cell_h, cell_ih, atoms(i,:), atoms(n_list(j),:), rij, r2) + IF (r2 < rc*rc) THEN ! Only calculates contributions between neighbouring particles. + sigmaij = 1 ! geometric average of interaction parameters + epsij = 1 + IF (i .le. n_type2) THEN + sigmaij = 0.6 + epsij = 0.6 + ENDIF + IF (i .le. n_type2) THEN + sigmaij = sigmaij*0.6 + epsij = epsij*0.6 + ENDIF + sigmaij = sigma*DSQRT(sigmaij) + epsij = eps*DSQRT(epsij) + + CALL LJ_fij(sigmaij, epsij, rij, sqrt(r2), pot_ij, fij) + + forces(i,:) = forces(i,:) + fij + forces(n_list(j),:) = forces(n_list(j),:) - fij + pot = pot + pot_ij + DO k = 1, 3 + DO l = k, 3 + ! Only the upper triangular elements calculated. + virial(k,l) = virial(k,l) + fij(k)*rij(l) + ENDDO + ENDDO + ENDIF + ENDDO + start = index_list(i) + 1 + ENDDO + + ! Assuming an upper-triangular vector matrix for the simulation box. + volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3) + CALL LJ_longrange(rc, sigma*(1-n_type2*0.4/natoms), eps*(1-n_type2*0.4/natoms), natoms, volume, pot_lr, vir_lr) + pot = pot + pot_lr + DO k = 1, 3 + virial(k,k) = virial(k,k) + vir_lr + ENDDO + + END SUBROUTINE END MODULE diff --git a/drivers/f90/driver.f90 b/drivers/f90/driver.f90 index 4ad8a4070..7777f69dd 100644 --- a/drivers/f90/driver.f90 +++ b/drivers/f90/driver.f90 @@ -175,6 +175,8 @@ PROGRAM DRIVER vstyle = 29 ELSEIF (trim(cmdbuffer) == "water_dip_pol") THEN vstyle = 31 + ELSEIF (trim(cmdbuffer) == "ljmix") THEN + vstyle = 32 ELSEIF (trim(cmdbuffer) == "qtip4pf-c-1") THEN vstyle = 60 ELSEIF (trim(cmdbuffer) == "qtip4pf-c-2") THEN @@ -193,7 +195,7 @@ PROGRAM DRIVER vstyle = 99 ! returns non-zero but otherwise meaningless values ELSE WRITE(*,*) " Unrecognized potential type ", trim(cmdbuffer) - WRITE(*,*) " Use -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4pf-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta] " + WRITE(*,*) " Use -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4pf-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta] " STOP "ENDED" ENDIF ELSEIF (ccmd == 4) THEN @@ -294,7 +296,6 @@ PROGRAM DRIVER STOP "ENDED" ENDIF isinit = .true. - ELSEIF (23 == vstyle) THEN !MB IF (par_count == 0) THEN ! defaults values vpars(1) = 0.004737803248674678 @@ -337,10 +338,10 @@ PROGRAM DRIVER STOP "ENDED" ENDIF isinit = .true. - ELSEIF (22 == vstyle) THEN !ljpolymer + ELSEIF (22 == vstyle .or. 32 == vstyle) THEN !ljpolymer or ljmix IF (4/= par_count) THEN WRITE(*,*) "Error: parameters not initialized correctly." - WRITE(*,*) "For ljpolymer potential use n_monomer,sigma,epsilon,cutoff" + WRITE(*,*) "For ljpolymer and ljmix potential use n_monomer,sigma,epsilon,cutoff" STOP "ENDED" ELSE n_monomer = nint(vpars(1)) @@ -901,6 +902,8 @@ PROGRAM DRIVER CALL SG_getall(rc, nat, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial) ELSEIF (vstyle == 22) THEN ! ljpolymer potential. CALL ljpolymer_getall(n_monomer,rc,sigma,eps,stiffness,nat,atoms,cell_h,cell_ih,index_list,n_list,pot,forces,virial) + ELSEIF (vstyle == 32) THEN ! lj mixture. + CALL ljmix_getall(n_monomer,rc,sigma,eps,nat,atoms,cell_h,cell_ih,index_list,n_list,pot,forces,virial) ENDIF IF (verbose > 0) WRITE(*,*) " Calculated energy is ", pot ENDIF @@ -1048,7 +1051,7 @@ PROGRAM DRIVER CONTAINS SUBROUTINE helpmessage ! Help banner - WRITE(*,*) " SYNTAX: driver.x [-u] -a address -p port -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4p-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta]" + WRITE(*,*) " SYNTAX: driver.x [-u] -a address -p port -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4p-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta]" WRITE(*,*) " -o 'comma_separated_parameters' [-v] " WRITE(*,*) "" WRITE(*,*) " For LJ potential use -o sigma,epsilon,cutoff " @@ -1056,7 +1059,7 @@ SUBROUTINE helpmessage WRITE(*,*) " For 1D/3D harmonic oscillator use -o k " WRITE(*,*) " For 1D morse oscillators use -o r0,D,a" WRITE(*,*) " For qtip4pf-efield use -o Ex,Ey,Ez with Ei in V/nm" - WRITE(*,*) " For ljpolymer use -o n_monomer,sigma,epsilon,cutoff " + WRITE(*,*) " For ljpolymer or lkmix use -o n_monomer,sigma,epsilon,cutoff " WRITE(*,*) " For gas, dummy, use the optional -o sleep_seconds to add a delay" WRITE(*,*) " For the ideal, qtip4pf*, zundel, ch4hcbe, nasa, doublewell or doublewell_1D no options are needed! " END SUBROUTINE helpmessage diff --git a/examples/features/monte_carlo/alchemical/init.xyz b/examples/features/monte_carlo/alchemical/init.xyz new file mode 100644 index 000000000..f8b5f9693 --- /dev/null +++ b/examples/features/monte_carlo/alchemical/init.xyz @@ -0,0 +1,8 @@ +6 +# CELL(abcABC): 130.18034 170.29206 240.21758 90.00000 90.00000 90.00000 cell{atomic_unit} Traj: positions{angstrom} Step: 12 Bead: 0 +H 1.3762859974402626 1.6475377505750781 1.3947309252535616 +C 0.5610492771673459 0.6910172235578630 0.5785825986233499 +H 0.0313332395759015 0.1474340996450444 1.3607366615876173 +D 1.3412729333952838 0.1418736998062487 0.0516893010969587 +D -0.0352236727959057 1.3786671745683925 -0.0209976311964991 +H -0.162629 -0.151689 -0.290930 diff --git a/examples/features/monte_carlo/alchemical/input.xml b/examples/features/monte_carlo/alchemical/input.xml new file mode 100644 index 000000000..f6f43e3a8 --- /dev/null +++ b/examples/features/monte_carlo/alchemical/input.xml @@ -0,0 +1,42 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, pressure_cv{megapascal} ] + positions{angstrom} + + + 1000 + + 31415 + + +
ch5p
+ + + + init.xyz + 100 + + + + + + + + 0.5 + + 10 + + + + + + 1 + [ H, D ] + + + + + 400 + + + diff --git a/examples/features/monte_carlo/alchemical/run.sh b/examples/features/monte_carlo/alchemical/run.sh new file mode 100644 index 000000000..853946533 --- /dev/null +++ b/examples/features/monte_carlo/alchemical/run.sh @@ -0,0 +1,16 @@ +ipi=i-pi +driver="i-pi-driver -m ch4hcbe -u -a ch5p" +sleep_time=4 + +${ipi} input.xml > log.i-pi & +echo "# i-PI is running" + +echo "# Waiting for ${sleep_time} (s) before executing driver" +sleep ${sleep_time} + +${driver} > /dev/null & +echo "# Driver is running" + +wait + +echo "# Simulation complete" diff --git a/examples/features/monte_carlo/atomswap/input.xml b/examples/features/monte_carlo/atomswap/input.xml new file mode 100644 index 000000000..c53c204eb --- /dev/null +++ b/examples/features/monte_carlo/atomswap/input.xml @@ -0,0 +1,42 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, pressure_cv{megapascal} ] + positions{angstrom} + + + 1000 + + 31415 + + +
ljmix
+ + + + ljmix.xyz + 100 + + + + + + + + 5 + + 10 + + + + + + 1 + [ Ne, Ar ] + + + + + 300 + + + diff --git a/examples/features/monte_carlo/atomswap/ljmix.xyz b/examples/features/monte_carlo/atomswap/ljmix.xyz new file mode 120000 index 000000000..d29475627 --- /dev/null +++ b/examples/features/monte_carlo/atomswap/ljmix.xyz @@ -0,0 +1 @@ +../../../init_files/ljmix.xyz \ No newline at end of file diff --git a/examples/features/monte_carlo/atomswap/run.sh b/examples/features/monte_carlo/atomswap/run.sh new file mode 100644 index 000000000..da35610eb --- /dev/null +++ b/examples/features/monte_carlo/atomswap/run.sh @@ -0,0 +1,16 @@ +ipi=i-pi +driver="i-pi-driver -m ljmix -u -a ljmix -o 32,4,1e-3,12" +sleep_time=4 + +${ipi} input.xml > log.i-pi & +echo "# i-PI is running" + +echo "# Waiting for ${sleep_time} (s) before executing driver" +sleep ${sleep_time} + +${driver} > /dev/null & +echo "# Driver is running" + +wait + +echo "# Simulation complete" diff --git a/examples/init_files/ljmix.xyz b/examples/init_files/ljmix.xyz new file mode 100644 index 000000000..f4e51001d --- /dev/null +++ b/examples/init_files/ljmix.xyz @@ -0,0 +1,110 @@ +108 +# CELL(abcABC): 17.8469984966 17.8469984966 17.8469984966 90.00000 90.00000 90.00000 Traj: positions{angstrom} Step: 0 Bead: 0 cell{angstrom} +Ne 0.147000137166 0.132000079971 0.114000117173 +Ne -0.125000123837 2.94299788633 3.18800164278 +Ne 3.04200165055 2.93700230854 0.107000161039 +Ne 2.97599737714 -0.147999752916 2.96000035008 +Ne -0.228999850096 -0.0300000085539 5.92297467609 +Ne 0.0939998645211 2.99900071046 8.91499553915 +Ne 2.95399948052 2.89500151338 5.91000983444 +Ne 2.8920010786 -0.155000238227 8.8700154763 +Ne -0.0480000242697 0.0869999083873 11.900983782 +Ne 0.0169999766244 2.94699846603 14.9550242141 +Ne 2.97899781192 2.98899926119 11.875001181 +Ne 2.97599737714 0.114000117173 14.7939955891 +Ne -0.0529997434675 5.85201201223 -0.165000099964 +Ne 0.00600002287786 8.85800315363 3.04000136069 +Ne 3.04100150562 8.95097958943 -0.0960001543749 +Ne 3.00700186988 5.87900004994 2.87599875977 +Ne 0.0569997939979 5.91000983444 6.03600692814 +Ne -0.0969997701246 8.8669991662 8.80397416049 +Ne 2.83199767476 8.76597923681 5.71897686163 +Ne 3.14800113748 5.95599533399 8.96399734879 +Ne 0.0309999947276 5.97398735913 11.9319935665 +Ne 0.12399997891 8.97198792467 14.7229800075 +Ne 2.99800056553 8.86101946373 11.9519964651 +Ne 2.97100194428 5.94297757463 14.9140129803 +Ne 0.00899998139907 11.7700124225 0.0890001982411 +Ne 0.0150000042769 14.9309995688 2.99199969597 +Ne 2.92400042449 14.8889828983 0.0619999894552 +Ne 2.83199767476 11.9180232882 3.16199787469 +Ne -0.0969997701246 11.8850026303 6.10501163633 +Ne 0.169000150494 14.7310235011 8.88901293814 +Ne 2.8769989047 14.7500209629 5.92101672041 +Ne 3.07800157614 12.0379877617 8.93600387439 +Ar 0.0190000018896 12.010999724 11.8850026303 +Ar 0.139000036105 14.8480245822 14.8409865253 +Ar 2.90399752595 14.8350068229 12.0199957365 +Ar 3.08300230078 12.0130105974 14.9420064547 +Ar 6.04198663062 0.050999982791 -0.122000218234 +Ar 5.91101527114 2.9569999153 2.81099992306 +Ar 8.94902163375 3.11100106696 0.0129999790117 +Ar 8.9909853865 -0.0100000204905 2.79399745931 +Ar 6.06098409246 0.0420000013919 6.018014903 +Ar 6.15602431937 3.00200114524 8.7079814146 +Ar 8.92198067832 3.0690002718 5.90900439774 +Ar 8.98802199413 0.132000079971 9.00601401927 +Ar 5.91598953692 0.156999998903 11.9360153133 +Ar 6.01097684611 2.8409989791 14.9089857968 +Ar 8.97399879806 2.98999940612 11.9060109655 +Ar 8.94600532365 -0.011999992838 14.7619803679 +Ar 5.87100947407 5.8559808413 -0.118999783453 +Ar 6.16401489524 8.84297452087 3.01699802737 +Ar 9.19302524528 8.76201040774 0.075999901723 +Ar 8.8339785083 5.96499134656 2.83499810954 +Ar 5.9870051185 6.00298627024 5.94101961895 +Ar 5.97700366923 8.97399879806 8.83900569179 +Ar 9.02501148111 8.84599083096 6.02198373208 +Ar 8.84800170436 5.69399969732 9.14301799894 +Ar 5.92398011279 6.03198518135 11.8649997317 +Ar 6.0190203397 8.82699336913 14.8660166074 +Ar 8.834983945 9.07602416415 11.9409895791 +Ar 8.89102381154 6.03902323824 14.8950155185 +Ar 5.8559808413 11.9979819646 0.104999871185 +Ar 5.82798736689 14.7829887031 2.90299738102 +Ar 8.94198357686 14.8989843475 -0.0420000013919 +Ar 8.93002417191 11.9660196611 3.01499773752 +Ar 6.07299641512 11.8649997317 5.91202070784 +Ar 6.03097974465 14.8909937717 8.97902598156 +Ar 8.97198792467 14.8580260315 5.88302179673 +Ar 8.95701220962 11.8349953839 8.96701365889 +Ar 5.96599678326 11.8819863202 12.0440203819 +Ar 5.94699932142 14.7049879824 14.7739926905 +Ar 8.87599517877 14.852998848 11.9509910284 +Ar 8.9919908232 11.9170178515 14.9289886954 +Ar 11.9460167626 0.0880000533142 -0.0279999832887 +Ar 11.8290156815 3.08100201092 3.11700193652 +Ar 15.0040260237 3.0330003462 -0.00799999522534 +Ar 14.8280216837 0.11999992838 3.0289997665 +Ar 11.9060109655 0.0169999766244 5.99499569437 +Ar 11.8729903076 2.91599926507 9.01601546854 +Ar 14.9069749234 3.06199925731 5.96001708079 +Ar 14.9479861572 -0.122000218234 8.8639828561 +Ar 11.8690214785 -0.118000167703 11.9299826931 +Ar 11.9319935665 2.86300216749 14.9789959417 +Ar 14.8350068229 2.97599737714 12.0179848631 +Ar 14.8519934113 -0.217999843432 14.7829887031 +Ar 11.8919877694 5.99002142859 0.0579999389248 +Ar 11.8159979221 8.68400968698 2.93600216361 +Ar 14.9579876065 8.87599517877 -0.0550000333213 +Ar 14.8399810886 6.06098409246 2.91599926507 +Ar 11.9109852313 5.97102396676 5.72601491853 +Ar 11.9209866806 8.94701076035 8.91499553915 +Ar 15.0320194981 8.97198792467 5.83100367699 +Ar 14.6990082799 5.88900149921 8.93097669089 +Ar 11.764985239 5.84402143636 11.8460022699 +Ar 11.9830062496 8.75502526857 14.9330104422 +Ar 14.8399810886 8.92600242512 11.9470221993 +Ar 14.9759796316 5.71199172246 14.7709763805 +Ar 11.823988498 11.8709794342 -0.0070000090516 +Ar 11.8939986428 14.8399810886 2.77099941776 +Ar 15.0390046373 14.9279832587 0.0209999742371 +Ar 14.8760180566 11.8509765357 2.90999839551 +Ar 11.7849881376 11.9340044399 6.07098554172 +Ar 11.93702075 14.6810162548 8.92901873521 +Ar 14.9350213156 14.9020006576 5.87302034746 +Ar 14.8819977591 11.9440058892 9.00701945597 +Ar 11.9409895791 11.9560182119 11.8770120544 +Ar 11.7910207578 14.8999897842 14.9550242141 +Ar 14.838023133 14.9469807205 11.8309736371 +Ar 14.8909937717 12.0050200215 15.0500115233 diff --git a/ipi/engine/motion/al6xxx_kmc.py b/ipi/engine/motion/al6xxx_kmc.py index b256574da..10ab09570 100755 --- a/ipi/engine/motion/al6xxx_kmc.py +++ b/ipi/engine/motion/al6xxx_kmc.py @@ -240,7 +240,7 @@ def bind(self, ens, beads, nm, cell, bforce, prng, omaker): idx = np.asarray( list(range(self.ncell**3)), int ) # initialize random distribution of atoms - self.prng.rng.shuffle(idx) + self.prng.shuffle(idx) self.idx = idx # initialize state based on the index diff --git a/ipi/engine/motion/alchemy.py b/ipi/engine/motion/alchemy.py index 995e80ae7..ee4448409 100644 --- a/ipi/engine/motion/alchemy.py +++ b/ipi/engine/motion/alchemy.py @@ -83,7 +83,7 @@ def AXlist(self, atomtype): def step(self, step=None): # picks number of attempted exchanges - ntries = self.prng.rng.poisson(self.nxc) + ntries = self.prng.poisson(self.nxc) if ntries == 0: return @@ -120,10 +120,10 @@ def step(self, step=None): # if (1.0/self.nxc < self.prng.u) : return # tries a round of exhanges with probability 1/nmc for x in range(ntries): - i = self.prng.rng.randint(lenlist) - j = self.prng.rng.randint(lenlist) + i = self.prng.integers(0, lenlist) + j = self.prng.integers(0, lenlist) while self.beads.names[axlist[i]] == self.beads.names[axlist[j]]: - j = self.prng.rng.randint(lenlist) # makes sure we pick a real exchange + j = self.prng.integers(0, lenlist) # makes sure we pick a real exchange # energy change due to the swap difspring = (atomspring[i] - atomspring[j]) * ( diff --git a/ipi/engine/motion/atomswap.py b/ipi/engine/motion/atomswap.py index e69c6978f..5ddcf0a54 100644 --- a/ipi/engine/motion/atomswap.py +++ b/ipi/engine/motion/atomswap.py @@ -81,7 +81,7 @@ def AXlist(self, atomtype): def step(self, step=None): # picks number of attempted exchanges - ntries = self.prng.rng.poisson(self.nxc) + ntries = self.prng.poisson(self.nxc) if ntries == 0: return @@ -104,10 +104,10 @@ def step(self, step=None): self.cell.h ) # just in case the cell gets updated in the other motion classes for x in range(ntries): - i = self.prng.rng.randint(lenlist) - j = self.prng.rng.randint(lenlist) + i = self.prng.integers(0, lenlist) + j = self.prng.integers(0, lenlist) while self.beads.names[axlist[i]] == self.beads.names[axlist[j]]: - j = self.prng.rng.randint(lenlist) # makes sure we pick a real exchange + j = self.prng.integers(0, lenlist) # makes sure we pick a real exchange # map the "subset" indices back to the "absolute" atom indices i = axlist[i] diff --git a/ipi/utils/prng.py b/ipi/utils/prng.py index b9f62c163..75c1fcf56 100644 --- a/ipi/utils/prng.py +++ b/ipi/utils/prng.py @@ -104,7 +104,7 @@ def g(self): return self.rng[0].standard_normal() - def gamma(self, k, theta=1.0): + def gamma(self, k, theta=1.0, size=None): """Interface to the standard gamma() function. Args: @@ -116,7 +116,43 @@ def gamma(self, k, theta=1.0): mean value theta. """ - return self.rng[0].gamma(k, theta) + return self.rng[0].gamma(k, theta, size) + + def poisson(self, lam=1.0, size=None): + """Interface to the standard poisson() function. + + Args: + lam: Mean of the Poisson distribution + + Returns: + A random number from a Poisson distribution + """ + + return self.rng[0].poisson(lam, size) + + def integers(self, low, high=None, size=None, dtype=np.int64, endpoint=False): + """Interface to the standard integers() function. + + Args: + Same as numpy.Generator.integers + + Returns: + Random integers in the prescribed interval + """ + + return self.rng[0].integers(low, high, size, dtype, endpoint) + + def shuffle(self, x, axis=0): + """Interface to the standard shuffle() function. + + Args: + Same as numpy.Generator.shuffle + + Returns: + None + """ + + self.rng[0].shuffle(x, axis) def gfill_serial(self, out): """Fills a pre-allocated array serially