From 45575394f98facbe12b9a92c0a81b75353eccb5d Mon Sep 17 00:00:00 2001 From: Bob Feng Date: Tue, 24 Oct 2017 14:06:13 -0700 Subject: [PATCH 1/2] added changes to super reptation in temporary file so old interpolation is not disturbed --- src/mc/MC_int_super_rep_temp.f03 | 337 +++++++++++++++++++++++++++++++ src/mc/MC_superReptation.f03 | 198 ++++++++++-------- 2 files changed, 452 insertions(+), 83 deletions(-) create mode 100644 src/mc/MC_int_super_rep_temp.f03 diff --git a/src/mc/MC_int_super_rep_temp.f03 b/src/mc/MC_int_super_rep_temp.f03 new file mode 100644 index 00000000..9909a759 --- /dev/null +++ b/src/mc/MC_int_super_rep_temp.f03 @@ -0,0 +1,337 @@ +!---------------------------------------------------------------! +! +! This subroutine calculates the change in the self energy for +! a super reptation move. That is a reptation move where the +! chain identities change along with position so that middle +! beads appear not to change. +! +! Created from MC_int_rep by Quinn on 8/10/17 +!--------------------------------------------------------------- +subroutine MC_int_super_rep_temp(wlc_p,wlc_d,IT1,IT2,forward) +use params,only: dp,wlcsim_params,wlcsim_data +implicit none + +! iputs +TYPE(wlcsim_params), intent(in) :: wlc_p ! <---- Contains output +TYPE(wlcsim_data), intent(inout) :: wlc_d +integer, intent(in) :: IT1 ! Test bead position 1 +integer, intent(in) :: IT2 ! Test bead position 2 + +! Internal variables +integer I ! For looping over bins +integer II ! For looping over IB +integer IB1,IB2 ! Bead index +integer rrdr ! -1 if r, 1 if r + dr +integer IX(2),IY(2),IZ(2) +real(dp) WX(2),WY(2),WZ(2) +real(dp) WTOT ! total weight ascribed to bin +real(dp) RBin(3) ! bead position +integer inDBin ! index of bin +integer ISX,ISY,ISZ +LOGICAL isA ! The bead is of type A +real(dp), dimension(-2:2) :: phi2 +integer m_index ! m from Ylm spherical harmonics +integer NBinX(3) +real(dp) temp !for speeding up code +LOGICAL, intent(in) :: forward ! move forward +integer AminusB +NBinX = wlc_p%NBinX + +wlc_d%NPHI = 0 + +! unroll the loop + +!(case #2: II = 2) +do IB2 = (IT2-wlc_p%nBPM+1),IT2 + if (forward) then + ! moving forward I2 is added + rrdr = 1 + else + ! moving forward I1 is removed + rrdr = -1 + endif + !else + ! print*, "Error in MC_int_rep, II = {1,2}" + ! stop 1 + !endif + ! subract current and add new + if (rrdr == -1) then + RBin(1) = wlc_d%R(1,IB2) + RBin(2) = wlc_d%R(2,IB2) + RBin(3) = wlc_d%R(3,IB2) + isA = wlc_d%AB(IB2).eq.1 + else + RBin(1) = wlc_d%RP(1,IB2) + RBin(2) = wlc_d%RP(2,IB2) + RBin(3) = wlc_d%RP(3,IB2) + isA = wlc_d%ABP(IB2).eq.1 + endif + if (wlc_p%chi_l2_on .and. isA) then + if (rrdr == -1) then + call Y2calc(wlc_d%U(:,IB2),phi2) + else + call Y2calc(wlc_d%UP(:,IB2),phi2) + endif + else + ! You could give some MS parameter to B as well if you wanted + phi2=0.0_dp + endif + ! -------------------------------------------------- + ! + ! Interpolate beads into bins + ! + ! -------------------------------------------------- + call interp(wlc_p,RBin,IX,IY,IZ,WX,WY,WZ) + + ! ------------------------------------------------------- + ! + ! Count beads in bins + ! + ! ------------------------------------------------------ + ! Add or Subtract volume fraction with weighting from each bin + ! I know that it looks bad to have this section of code twice but it + ! makes it faster. + if (isA) then + do ISX = 1,2 + if ((IX(ISX).le.0).OR.(IX(ISX).ge.(NBinX(1) + 1))) CYCLE + do ISY = 1,2 + if ((IY(ISY).le.0).OR.(IY(ISY).ge.(NBinX(2) + 1))) CYCLE + do ISZ = 1,2 + if ((IZ(ISZ).le.0).OR.(IZ(ISZ).ge.(NBinX(3) + 1))) cycle + WTOT = WX(ISX)*WY(ISY)*WZ(ISZ) + inDBin = IX(ISX) + (IY(ISY)-1)*NBinX(1) + (IZ(ISZ)-1)*NBinX(1)*NBinX(2) + ! Generate list of which phi's change and by how much + I = wlc_d%NPHI + do + if (I.eq.0) then + wlc_d%NPHI = wlc_d%NPHI + 1 + wlc_d%inDPHI(wlc_d%NPHI) = inDBin + temp = rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + wlc_d%DPHIA(wlc_d%NPHI) = temp + wlc_d%DPHIB(wlc_d%NPHI) = 0.0_dp + if(wlc_p%chi_l2_on) then + do m_index = -2,2 + wlc_d%DPHI_l2(m_index,wlc_d%NPHI) = & + + phi2(m_index)*temp + enddo + endif + exit + elseif (inDBin == wlc_d%inDPHI(I)) then + temp = rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + wlc_d%DPHIA(I) = wlc_d%DPHIA(I) + temp + if(wlc_p%chi_l2_on) then + do m_index = -2,2 + wlc_d%DPHI_l2(m_index,I) = wlc_d%DPHI_l2(m_index,I) & + + phi2(m_index)*temp + enddo + endif + exit + else + I = I-1 + endif + enddo + enddo + enddo + enddo + else ! (isB) + do ISX = 1,2 + if ((IX(ISX).le.0).OR.(IX(ISX).ge.(NBinX(1) + 1))) CYCLE + do ISY = 1,2 + if ((IY(ISY).le.0).OR.(IY(ISY).ge.(NBinX(2) + 1))) CYCLE + do ISZ = 1,2 + if ((IZ(ISZ).le.0).OR.(IZ(ISZ).ge.(NBinX(3) + 1))) cycle + WTOT = WX(ISX)*WY(ISY)*WZ(ISZ) + inDBin = IX(ISX) + (IY(ISY)-1)*NBinX(1) + (IZ(ISZ)-1)*NBinX(1)*NBinX(2) + ! Generate list of which phi's change and by how much + I = wlc_d%NPHI + do + if (I.eq.0) then + wlc_d%NPHI = wlc_d%NPHI + 1 + wlc_d%inDPHI(wlc_d%NPHI) = inDBin + wlc_d%DPHIA(wlc_d%NPHI) = 0.0_dp + wlc_d%DPHIB(wlc_d%NPHI) = rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + if(wlc_p%chi_l2_on) then + do m_index = -2,2 + ! This is somewhat wastefull, could eliminate for speedup by having another NPHI for L=2 + wlc_d%DPHI_l2(m_index,wlc_d%NPHI) = 0.0_dp + enddo + endif + exit + elseif (inDBin == wlc_d%inDPHI(I)) then + wlc_d%DPHIB(I) = wlc_d%DPHIB(I) + rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + exit + else + I = I-1 + endif + enddo + enddo !ISZ + enddo !ISY + enddo !ISX + endif +enddo ! loop over IB A.k.a. beads +!---------------------------------------------------------------! + +!(case # 1: II = 1) + +do IB1 = IT1,(IT1 + wlc_p%nBPM - 1) + ! if (II.eq.1) then + ! IB = I1 + if (forward) then + ! moving forward I1 is removed + rrdr = -1 + else + ! moving backward I2 is added + rrdr = 1 + endif +! elseif (II.eq.2) then + ! IB = I2 + ! if (forward) then + ! ! moving forward I2 is added + ! rrdr = 1 + ! else + ! moving forward I1 is removed + ! rrdr = -1 + ! endif + !else + ! print*, "Error in MC_int_rep, II = {1,2}" + ! stop 1 + !endif + ! subract current and add new + if (rrdr == -1) then + RBin(1) = wlc_d%R(1,IB1) + RBin(2) = wlc_d%R(2,IB1) + RBin(3) = wlc_d%R(3,IB1) + isA = wlc_d%AB(IB1).eq.1 + else + RBin(1) = wlc_d%RP(1,IB1) + RBin(2) = wlc_d%RP(2,IB1) + RBin(3) = wlc_d%RP(3,IB1) + isA = wlc_d%ABP(IB1).eq.1 + endif + if (wlc_p%chi_l2_on .and. isA) then + if (rrdr == -1) then + call Y2calc(wlc_d%U(:,IB1),phi2) + else + call Y2calc(wlc_d%UP(:,IB1),phi2) + endif + else + ! You could give some MS parameter to B as well if you wanted + phi2=0.0_dp + endif + ! -------------------------------------------------- + ! + ! Interpolate beads into bins + ! + ! -------------------------------------------------- + call interp(wlc_p,RBin,IX,IY,IZ,WX,WY,WZ) + + ! ------------------------------------------------------- + ! + ! Count beads in bins + ! + ! ------------------------------------------------------ + ! Add or Subtract volume fraction with weighting from each bin + ! I know that it looks bad to have this section of code twice but it + ! makes it faster. + if (isA) then + do ISX = 1,2 + if ((IX(ISX).le.0).OR.(IX(ISX).ge.(NBinX(1) + 1))) CYCLE + do ISY = 1,2 + if ((IY(ISY).le.0).OR.(IY(ISY).ge.(NBinX(2) + 1))) CYCLE + do ISZ = 1,2 + if ((IZ(ISZ).le.0).OR.(IZ(ISZ).ge.(NBinX(3) + 1))) cycle + WTOT = WX(ISX)*WY(ISY)*WZ(ISZ) + inDBin = IX(ISX) + (IY(ISY)-1)*NBinX(1) + (IZ(ISZ)-1)*NBinX(1)*NBinX(2) + ! Generate list of which phi's change and by how much + I = wlc_d%NPHI + do + if (I.eq.0) then + wlc_d%NPHI = wlc_d%NPHI + 1 + wlc_d%inDPHI(wlc_d%NPHI) = inDBin + temp = rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + wlc_d%DPHIA(wlc_d%NPHI) = temp + wlc_d%DPHIB(wlc_d%NPHI) = 0.0_dp + if(wlc_p%chi_l2_on) then + do m_index = -2,2 + wlc_d%DPHI_l2(m_index,wlc_d%NPHI) = & + + phi2(m_index)*temp + enddo + endif + exit + elseif (inDBin == wlc_d%inDPHI(I)) then + temp = rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + wlc_d%DPHIA(I) = wlc_d%DPHIA(I) + temp + if(wlc_p%chi_l2_on) then + do m_index = -2,2 + wlc_d%DPHI_l2(m_index,I) = wlc_d%DPHI_l2(m_index,I) & + + phi2(m_index)*temp + enddo + endif + exit + else + I = I-1 + endif + enddo + enddo + enddo + enddo + else ! isB + do ISX = 1,2 + if ((IX(ISX).le.0).OR.(IX(ISX).ge.(NBinX(1) + 1))) CYCLE + do ISY = 1,2 + if ((IY(ISY).le.0).OR.(IY(ISY).ge.(NBinX(2) + 1))) CYCLE + do ISZ = 1,2 + if ((IZ(ISZ).le.0).OR.(IZ(ISZ).ge.(NBinX(3) + 1))) cycle + WTOT = WX(ISX)*WY(ISY)*WZ(ISZ) + inDBin = IX(ISX) + (IY(ISY)-1)*NBinX(1) + (IZ(ISZ)-1)*NBinX(1)*NBinX(2) + ! Generate list of which phi's change and by how much + I = wlc_d%NPHI + do + if (I.eq.0) then + wlc_d%NPHI = wlc_d%NPHI + 1 + wlc_d%inDPHI(wlc_d%NPHI) = inDBin + wlc_d%DPHIA(wlc_d%NPHI) = 0.0_dp + wlc_d%DPHIB(wlc_d%NPHI) = rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + if(wlc_p%chi_l2_on) then + do m_index = -2,2 + ! This is somewhat wastefull, could eliminate for speedup by having another NPHI for L=2 + wlc_d%DPHI_l2(m_index,wlc_d%NPHI) = 0.0_dp + enddo + endif + exit + elseif (inDBin == wlc_d%inDPHI(I)) then + wlc_d%DPHIB(I) = wlc_d%DPHIB(I) + rrdr*WTOT*wlc_p%beadVolume/wlc_d%Vol(inDBin) + exit + else + I = I-1 + endif + enddo + enddo !ISZ + enddo !ISY + enddo !ISX + endif +enddo! loop over IB A.k.a. beads +! --------------------------------------------------------------------- +! + + + +!---------------------------------------------------------- +! +! Intermediate Beads Don't change field +! +!----------------------------------------------------------- + +! ... skipping + +! --------------------------------------------------------------------- +! +! Calcualte change in energy +! +!--------------------------------------------------------------------- +call hamiltonian(wlc_p,wlc_d,.false.) + +!RETURN +END subroutine + +!---------------------------------------------------------------! diff --git a/src/mc/MC_superReptation.f03 b/src/mc/MC_superReptation.f03 index 31bfc996..ddb1c4a0 100644 --- a/src/mc/MC_superReptation.f03 +++ b/src/mc/MC_superReptation.f03 @@ -40,6 +40,9 @@ subroutine MC_superReptation(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& real(dp) r_relative(3) ! r in new coordinate system real(dp) u_relative(3) ! u in new coordinate system logical, intent(out) :: forward +integer beadi +real(dp), allocatable :: Rtemp(:,:) +real(dp), allocatable :: Utemp(:,:) !TOdo saving RP is not actually needed, even in these cases, but Brad's code assumes that we have RP. if (wlc_p%ring .OR. wlc_p%interp_bead_lennard_jones) then @@ -47,70 +50,86 @@ subroutine MC_superReptation(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& UP = U endif +allocate(Rtemp(3,wlc_p%nb)) +allocate(Utemp(3,wlc_p%nb)) + + ! single bead reptation call random_index(wlc_p%NP,irnd,rand_stat) IP=irnd(1) IT1 = wlc_p%NB*(IP-1) + 1 IT2 = wlc_p%NB*(IP-1) + wlc_p%NB -IB1 = mod(IT1,wlc_p%NB) -IB2 = mod(IT2,wlc_p%NB) + + +IB1 = 1 +IB2 = wlc_p%NB + +Rtemp(1:3,IB1:IB2) = R(1:3,IT1:IT2) +Utemp(1:3,IB1:IB2) = U(1:3,IT1:IT2) + ! move forward or backward call random_number(urnd,rand_stat) if (urnd(1).lt.0.5_dp) then forward = .true. - dR(1) = R(1,IT1 + 1)-R(1,IT1) - dR(2) = R(2,IT1 + 1)-R(2,IT1) - dR(3) = R(3,IT1 + 1)-R(3,IT1) - - Uvec(1) = U(1,IT1); Uvec(2) = U(2,IT1); Uvec(3) = U(3,IT1) - ! chose coordinate system - call random_perp(Uvec,pDir,tDir,rand_stat) - ! find next r and u in new coordinate system - u_relative(1) = Uvec(1)*U(1,IT1 + 1) + & - Uvec(2)*U(2,IT1 + 1) + & - Uvec(3)*U(3,IT1 + 1) - u_relative(2) = pDir(1)*U(1,IT1 + 1) + & - pDir(2)*U(2,IT1 + 1) + & - pDir(3)*U(3,IT1 + 1) - u_relative(3) = tDir(1)*U(1,IT1 + 1) + & - tDir(2)*U(2,IT1 + 1) + & - tDir(3)*U(3,IT1 + 1) - r_relative(1) = Uvec(1)*dR(1) + & - Uvec(2)*dR(2) + & - Uvec(3)*dR(3) - r_relative(2) = pDir(1)*dR(1) + & - pDir(2)*dR(2) + & - pDir(3)*dR(3) - r_relative(3) = tDir(1)*dR(1) + & - tDir(2)*dR(2) + & - tDir(3)*dR(3) - - - ! orient coordinate system with end of chain - Uvec(1) = U(1,IT2); Uvec(2) = U(2,IT2); Uvec(3) = U(3,IT2) - call random_perp(Uvec,pDir,tDir,rand_stat) - ! update UP and RP - UP(1,IT2) = Uvec(1)*u_relative(1) + pDir(1)*u_relative(2) + tDir(1)*u_relative(3) - UP(2,IT2) = Uvec(2)*u_relative(1) + pDir(2)*u_relative(2) + tDir(2)*u_relative(3) - UP(3,IT2) = Uvec(3)*u_relative(1) + pDir(3)*u_relative(2) + tDir(3)*u_relative(3) - mag = sqrt(UP(1,IT2)**2 + UP(2,IT2)**2 + UP(3,IT2)**2) - UP(1,IT2) = UP(1,IT2)/mag - UP(2,IT2) = UP(2,IT2)/mag - UP(3,IT2) = UP(3,IT2)/mag - RP(1,IT2) = R(1,IT2) + Uvec(1)*r_relative(1) + pDir(1)*r_relative(2) + tDir(1)*r_relative(3) - RP(2,IT2) = R(2,IT2) + Uvec(2)*r_relative(1) + pDir(2)*r_relative(2) + tDir(2)*r_relative(3) - RP(3,IT2) = R(3,IT2) + Uvec(3)*r_relative(1) + pDir(3)*r_relative(2) + tDir(3)*r_relative(3) - - do I = IT1,IT2-1 - RP(1,I) = R(1,I + 1) - RP(2,I) = R(2,I + 1) - RP(3,I) = R(3,I + 1) - UP(1,I) = U(1,I + 1) - UP(2,I) = U(2,I + 1) - UP(3,I) = U(3,I + 1) - ABP(I) = AB(I+1) + ! note, this is not the most efficient way of doing things + ! it would be more efficient to move all nbpm at a time. + do beadi =1,wlc_p%nbpm + dR(1) = Rtemp(1,IB1 + 1)-Rtemp(1,IB1) + dR(2) = Rtemp(2,IB1 + 1)-Rtemp(2,IB1) + dR(3) = Rtemp(3,IB1 + 1)-Rtemp(3,IB1) + + Uvec(1) = Utemp(1,IB1) + Uvec(2) = Utemp(2,IB1) + Uvec(3) = Utemp(3,IB1) + ! chose coordinate system + call random_perp(Uvec,pDir,tDir,rand_stat) + ! find next r and u in new coordinate system + u_relative(1) = Uvec(1)*Utemp(1,IB1 + 1) + & + Uvec(2)*Utemp(2,IB1 + 1) + & + Uvec(3)*Utemp(3,IB1 + 1) + u_relative(2) = pDir(1)*Utemp(1,IB1 + 1) + & + pDir(2)*Utemp(2,IB1 + 1) + & + pDir(3)*Utemp(3,IB1 + 1) + u_relative(3) = tDir(1)*Utemp(1,IB1 + 1) + & + tDir(2)*Utemp(2,IB1 + 1) + & + tDir(3)*Utemp(3,IB1 + 1) + r_relative(1) = Uvec(1)*dR(1) + & + Uvec(2)*dR(2) + & + Uvec(3)*dR(3) + r_relative(2) = pDir(1)*dR(1) + & + pDir(2)*dR(2) + & + pDir(3)*dR(3) + r_relative(3) = tDir(1)*dR(1) + & + tDir(2)*dR(2) + & + tDir(3)*dR(3) + + + ! orient coordinate system with end of chain + Uvec(1) = Utemp(1,IB2) + Uvec(2) = Utemp(2,IB2) + Uvec(3) = Utemp(3,IB2) + call random_perp(Uvec,pDir,tDir,rand_stat) + ! update UP and RP + UP(1,IT2) = Uvec(1)*u_relative(1) + pDir(1)*u_relative(2) + tDir(1)*u_relative(3) + UP(2,IT2) = Uvec(2)*u_relative(1) + pDir(2)*u_relative(2) + tDir(2)*u_relative(3) + UP(3,IT2) = Uvec(3)*u_relative(1) + pDir(3)*u_relative(2) + tDir(3)*u_relative(3) + mag = sqrt(UP(1,IT2)**2 + UP(2,IT2)**2 + UP(3,IT2)**2) + UP(1,IT2) = UP(1,IT2)/mag + UP(2,IT2) = UP(2,IT2)/mag + UP(3,IT2) = UP(3,IT2)/mag + RP(1,IT2) = Rtemp(1,IB2) + Uvec(1)*r_relative(1) + pDir(1)*r_relative(2) + tDir(1)*r_relative(3) + RP(2,IT2) = Rtemp(2,IB2) + Uvec(2)*r_relative(1) + pDir(2)*r_relative(2) + tDir(2)*r_relative(3) + RP(3,IT2) = Rtemp(3,IB2) + Uvec(3)*r_relative(1) + pDir(3)*r_relative(2) + tDir(3)*r_relative(3) + + RP(1:3,IT1:(IT2-1))=Rtemp(1:3,2:IB2) + UP(1:3,IT1:(IT2-1))=Utemp(1:3,2:IB2) + + ! pointing Rtemp to RP saves time compared to defining Rtemp as RP every loop. + Rtemp(1:3,IB1:IB2) = RP(1:3,IT1:IT2) + Utemp(1:3,IB1:IB2) = UP(1:3,IT1:IT2) enddo - ABP(IT2)=AB(IT1) ! extend cheical sequence + ABP(IT1:(IT2-wlc_p%nbpm))=AB((IT1+wlc_p%nbpm):IT2) + ABP((IT2-wlc_p%nbpm+1):IT2)=AB(IT1:(IT1+wlc_p%nbpm-1)) ! put end segment type on other end for detail balance ! RperpMag = sqrt(r_relative(2)**2 + r_relative(3)**2) ! RparaMag = r_relative(1) @@ -118,24 +137,25 @@ subroutine MC_superReptation(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& else forward = .false. - dR(1) = R(1,IT2)-R(1,IT2-1) - dR(2) = R(2,IT2)-R(2,IT2-1) - dR(3) = R(3,IT2)-R(3,IT2-1) + do beadi = 1,wlc_p%nbpm + dR(1) = Rtemp(1,IB2)-Rtemp(1,IB2-1) + dR(2) = Rtemp(2,IB2)-Rtemp(2,IB2-1) + dR(3) = Rtemp(3,IB2)-Rtemp(3,IB2-1) - Uvec(1) = U(1,IT2); Uvec(2) = U(2,IT2); Uvec(3) = U(3,IT2) + Uvec(1) = Utemp(1,IB2); Uvec(2) = Utemp(2,IB2); Uvec(3) = Utemp(3,IB2) ! chose coordinate system call random_perp(Uvec,pDir,tDir,rand_stat) ! find next r and u in new coordinate system - u_relative(1) = Uvec(1)*U(1,IT2-1) + & - Uvec(2)*U(2,IT2-1) + & - Uvec(3)*U(3,IT2-1) - u_relative(2) = pDir(1)*U(1,IT2-1) + & - pDir(2)*U(2,IT2-1) + & - pDir(3)*U(3,IT2-1) - u_relative(3) = tDir(1)*U(1,IT2-1) + & - tDir(2)*U(2,IT2-1) + & - tDir(3)*U(3,IT2-1) + u_relative(1) = Uvec(1)*Utemp(1,IB2-1) + & + Uvec(2)*Utemp(2,IB2-1) + & + Uvec(3)*Utemp(3,IB2-1) + u_relative(2) = pDir(1)*Utemp(1,IB2-1) + & + pDir(2)*Utemp(2,IB2-1) + & + pDir(3)*Utemp(3,IB2-1) + u_relative(3) = tDir(1)*Utemp(1,IB2-1) + & + tDir(2)*Utemp(2,IB2-1) + & + tDir(3)*Utemp(3,IB2-1) r_relative(1) = Uvec(1)*dR(1) + & Uvec(2)*dR(2) + & Uvec(3)*dR(3) @@ -147,7 +167,7 @@ subroutine MC_superReptation(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& tDir(3)*dR(3) ! orient coordinate system with end of chain - Uvec(1) = U(1,IT1); Uvec(2) = U(2,IT1); Uvec(3) = U(3,IT1) + Uvec(1) = Utemp(1,IB1); Uvec(2) = Utemp(2,IB1); Uvec(3) = Utemp(3,IB1) call random_perp(Uvec,pDir,tDir,rand_stat) ! update UP and RP UP(1,IT1) = Uvec(1)*u_relative(1) + pDir(1)*u_relative(2) + tDir(1)*u_relative(3) @@ -157,20 +177,32 @@ subroutine MC_superReptation(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& UP(1,IT1) = UP(1,IT1)/mag UP(2,IT1) = UP(2,IT1)/mag UP(3,IT1) = UP(3,IT1)/mag - RP(1,IT1) = R(1,IT1)-Uvec(1)*r_relative(1)-pDir(1)*r_relative(2)-tDir(1)*r_relative(3) - RP(2,IT1) = R(2,IT1)-Uvec(2)*r_relative(1)-pDir(2)*r_relative(2)-tDir(2)*r_relative(3) - RP(3,IT1) = R(3,IT1)-Uvec(3)*r_relative(1)-pDir(3)*r_relative(2)-tDir(3)*r_relative(3) - - do I = IT1 + 1,IT2 - RP(1,I) = R(1,I-1) - RP(2,I) = R(2,I-1) - RP(3,I) = R(3,I-1) - UP(1,I) = U(1,I-1) - UP(2,I) = U(2,I-1) - UP(3,I) = U(3,I-1) - ABP(I) = AB(I-1) - enddo - ABP(IT1)=AB(IT1) ! extend chemical sequence -endif + RP(1,IT1) = Rtemp(1,IB1)-Uvec(1)*r_relative(1)-pDir(1)*r_relative(2)-tDir(1)*r_relative(3) + RP(2,IT1) = Rtemp(2,IB1)-Uvec(2)*r_relative(1)-pDir(2)*r_relative(2)-tDir(2)*r_relative(3) + RP(3,IT1) = Rtemp(3,IB1)-Uvec(3)*r_relative(1)-pDir(3)*r_relative(2)-tDir(3)*r_relative(3) + + + RP(1:3,(IT1+1):IT2) = Rtemp(1:3,1:(IB2-1)) + UP(1:3,(IT1+1):IT2) = Utemp(1:3,1:(IB2-1)) + + Rtemp(1:3,IB1:IB2) = RP(1:3,IT1:IT2) + Utemp(1:3,IB1:IB2) = UP(1:3,IT1:IT2) + enddo + !do I = IT1 + 1,IT2 + ! RP(1,I) = R(1,I-1) + ! RP(2,I) = R(2,I-1) + ! RP(3,I) = R(3,I-1) + !UP(1,I) = U(1,I-1) + ! UP(2,I) = U(2,I-1) + ! UP(3,I) = U(3,I-1) + ! ABP(I) = AB(I-1) + ! enddo + !ABP(IT1)=AB(IT1) ! extend chemical sequence + + ABP((IT1+wlc_p%nbpm):IT2)=AB(IT1:(IT2-wlc_p%nbpm)) + ABP(IT1:(IT1+wlc_p%nbpm-1))=AB((IT2-wlc_p%nbpm+1):IT2) ! put end segment type on other end for detail balance + endif + deallocate(Rtemp) + deallocate(Utemp) end subroutine From 0c0491b944b37ee70ff8c2c6e9caf0519c8845f5 Mon Sep 17 00:00:00 2001 From: Bob Feng Date: Tue, 24 Oct 2017 14:15:38 -0700 Subject: [PATCH 2/2] super reptation in temporary file --- src/mc/MC_superReptationMulti.f03 | 103 +++++++++++++++++------------- 1 file changed, 58 insertions(+), 45 deletions(-) diff --git a/src/mc/MC_superReptationMulti.f03 b/src/mc/MC_superReptationMulti.f03 index ac637d8d..a7f6788e 100644 --- a/src/mc/MC_superReptationMulti.f03 +++ b/src/mc/MC_superReptationMulti.f03 @@ -15,12 +15,12 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& implicit none type(wlcsim_params), intent(in) :: wlc_p -real(dp), intent(in) :: R(3,wlc_p%NT) ! Bead positions -real(dp), intent(in) :: U(3,wlc_p%NT) ! Tangent vectors +real(dp), intent(in),target :: R(3,wlc_p%NT) ! Bead positions +real(dp), intent(in),target :: U(3,wlc_p%NT) ! Tangent vectors integer, intent(in) :: AB(wlc_p%NT) ! Current chemical identity integer, intent(out) :: ABP(wlc_p%NT) ! Proposed chemical identity -real(dp), intent(out) :: RP(3,wlc_p%NT) ! Bead positions -real(dp), intent(out) :: UP(3,wlc_p%NT) ! Tangent vectors +real(dp), intent(out),target :: RP(3,wlc_p%NT) ! Bead positions +real(dp), intent(out),target :: UP(3,wlc_p%NT) ! Tangent vectors integer, intent(out) :: IP ! Test polymer integer, intent(out) :: IT1 ! Index of test bead 1 integer, intent(out) :: IT2 ! Index of test bead 2 @@ -41,8 +41,8 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& real(dp) u_relative(3) ! u in new coordinate system logical, intent(out) :: forward integer beadi -real(dp), allocatable :: Rtemp(:,:) -real(dp), allocatable :: Utemp(:,:) +real(dp), allocatable :: Rtemp(:,:) +real(dp), allocatable :: Utemp(:,:) !TOdo saving RP is not actually needed, even in these cases, but Brad's code assumes that we have RP. if (wlc_p%ring .OR. wlc_p%interp_bead_lennard_jones) then @@ -50,8 +50,8 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& UP = U endif -allocate( Rtemp(3,wlc_p%nb) ) -allocate( Utemp(3,wlc_p%nb) ) +allocate(Rtemp(3,wlc_p%nb)) +allocate(Utemp(3,wlc_p%nb)) ! single bead reptation @@ -60,11 +60,13 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& IT1 = wlc_p%NB*(IP-1) + 1 IT2 = wlc_p%NB*(IP-1) + wlc_p%NB -Rtemp=R(:,IT1:IT2) -Utemp=U(:,IT1:IT2) IB1 = 1 IB2 = wlc_p%NB + +Rtemp(1:3,IB1:IB2) = R(1:3,IT1:IT2) +Utemp(1:3,IB1:IB2) = U(1:3,IT1:IT2) + ! move forward or backward call random_number(urnd,rand_stat) if (urnd(1).lt.0.5_dp) then @@ -76,8 +78,8 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& dR(2) = Rtemp(2,IB1 + 1)-Rtemp(2,IB1) dR(3) = Rtemp(3,IB1 + 1)-Rtemp(3,IB1) - Uvec(1) = Utemp(1,IB1); - Uvec(2) = Utemp(2,IB1); + Uvec(1) = Utemp(1,IB1) + Uvec(2) = Utemp(2,IB1) Uvec(3) = Utemp(3,IB1) ! chose coordinate system call random_perp(Uvec,pDir,tDir,rand_stat) @@ -103,8 +105,8 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& ! orient coordinate system with end of chain - Uvec(1) = Utemp(1,IB2); - Uvec(2) = Utemp(2,IB2); + Uvec(1) = Utemp(1,IB2) + Uvec(2) = Utemp(2,IB2) Uvec(3) = Utemp(3,IB2) call random_perp(Uvec,pDir,tDir,rand_stat) ! update UP and RP @@ -120,8 +122,10 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& RP(3,IT2) = Rtemp(3,IB2) + Uvec(3)*r_relative(1) + pDir(3)*r_relative(2) + tDir(3)*r_relative(3) RP(:,IT1:IT2-1)=Rtemp(:,2:IB2) - Rtemp=RP(:,IT1:IT2) - Utemp=UP(:,IT1:IT2) + + ! pointing Rtemp to RP saves time compared to defining Rtemp as RP every loop. + Rtemp(:,IB1:IB2) = RP(:,IT1:IT2) + Utemp(:,IB1:IB2) = UP(:,IT1:IT2) enddo ABP(IT1:IT2-wlc_p%nbpm)=ABP(IT1+wlc_p%nbpm:IT2) ABP(IT2-wlc_p%nbpm+1:IT2)=AB(IT1:IT1+wlc_p%nbpm-1) ! put end segment type on other end for detail balance @@ -132,24 +136,25 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& else forward = .false. - dR(1) = R(1,IT2)-R(1,IT2-1) - dR(2) = R(2,IT2)-R(2,IT2-1) - dR(3) = R(3,IT2)-R(3,IT2-1) + do beadi = 1,wlc_p%nbpm + dR(1) = Rtemp(1,IB2)-R(1,IB2-1) + dR(2) = Rtemp(2,IB2)-R(2,IB2-1) + dR(3) = Rtemp(3,IB2)-R(3,IB2-1) - Uvec(1) = U(1,IT2); Uvec(2) = U(2,IT2); Uvec(3) = U(3,IT2) + Uvec(1) = Utemp(1,IB2); Uvec(2) = Utemp(2,IB2); Uvec(3) = Utemp(3,IB2) ! chose coordinate system call random_perp(Uvec,pDir,tDir,rand_stat) ! find next r and u in new coordinate system - u_relative(1) = Uvec(1)*U(1,IT2-1) + & - Uvec(2)*U(2,IT2-1) + & - Uvec(3)*U(3,IT2-1) - u_relative(2) = pDir(1)*U(1,IT2-1) + & - pDir(2)*U(2,IT2-1) + & - pDir(3)*U(3,IT2-1) - u_relative(3) = tDir(1)*U(1,IT2-1) + & - tDir(2)*U(2,IT2-1) + & - tDir(3)*U(3,IT2-1) + u_relative(1) = Uvec(1)*Utemp(1,IB2-1) + & + Uvec(2)*Utemp(2,IB2-1) + & + Uvec(3)*Utemp(3,IB2-1) + u_relative(2) = pDir(1)*Utemp(1,IB2-1) + & + pDir(2)*Utemp(2,IB2-1) + & + pDir(3)*Utemp(3,IB2-1) + u_relative(3) = tDir(1)*Utemp(1,IB2-1) + & + tDir(2)*Utemp(2,IB2-1) + & + tDir(3)*Utemp(3,IB2-1) r_relative(1) = Uvec(1)*dR(1) + & Uvec(2)*dR(2) + & Uvec(3)*dR(3) @@ -161,7 +166,7 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& tDir(3)*dR(3) ! orient coordinate system with end of chain - Uvec(1) = U(1,IT1); Uvec(2) = U(2,IT1); Uvec(3) = U(3,IT1) + Uvec(1) = U(1,IB1); Uvec(2) = U(2,IB1); Uvec(3) = U(3,IB1) call random_perp(Uvec,pDir,tDir,rand_stat) ! update UP and RP UP(1,IT1) = Uvec(1)*u_relative(1) + pDir(1)*u_relative(2) + tDir(1)*u_relative(3) @@ -171,22 +176,30 @@ subroutine MC_superReptationMulti(wlc_p,R,U,RP,UP,AB,ABP,IP,IT1,IT2,IB1,IB2& UP(1,IT1) = UP(1,IT1)/mag UP(2,IT1) = UP(2,IT1)/mag UP(3,IT1) = UP(3,IT1)/mag - RP(1,IT1) = R(1,IT1)-Uvec(1)*r_relative(1)-pDir(1)*r_relative(2)-tDir(1)*r_relative(3) - RP(2,IT1) = R(2,IT1)-Uvec(2)*r_relative(1)-pDir(2)*r_relative(2)-tDir(2)*r_relative(3) - RP(3,IT1) = R(3,IT1)-Uvec(3)*r_relative(1)-pDir(3)*r_relative(2)-tDir(3)*r_relative(3) - - do I = IT1 + 1,IT2 - RP(1,I) = R(1,I-1) - RP(2,I) = R(2,I-1) - RP(3,I) = R(3,I-1) - UP(1,I) = U(1,I-1) - UP(2,I) = U(2,I-1) - UP(3,I) = U(3,I-1) - ABP(I) = AB(I-1) + RP(1,IT1) = R(1,IB1)-Uvec(1)*r_relative(1)-pDir(1)*r_relative(2)-tDir(1)*r_relative(3) + RP(2,IT1) = R(2,IB1)-Uvec(2)*r_relative(1)-pDir(2)*r_relative(2)-tDir(2)*r_relative(3) + RP(3,IT1) = R(3,IB1)-Uvec(3)*r_relative(1)-pDir(3)*r_relative(2)-tDir(3)*r_relative(3) + + + + RP(:,IT1+1:IT2) = Rtemp(:,1:IB2-1) + Rtemp(:,IB1:IB2) = RP(:,IT1:IT2) + Utemp(:,IB1:IB2) = UP(:,IT1:IT2) enddo - ABP(IT1)=AB(IT1) ! extend chemical sequence + !do I = IT1 + 1,IT2 + ! RP(1,I) = R(1,I-1) + ! RP(2,I) = R(2,I-1) + ! RP(3,I) = R(3,I-1) + !UP(1,I) = U(1,I-1) + ! UP(2,I) = U(2,I-1) + ! UP(3,I) = U(3,I-1) + ! ABP(I) = AB(I-1) + ! enddo + !ABP(IT1)=AB(IT1) ! extend chemical sequence endif -deallocate(Rtemp) -deallocate(Utemp) + ABP(IT1+wlc_p%nbpm:IT2)=ABP(IT1:IT2-wlc_p%nbpm) + ABP(IT1:IT1+wlc_p%nbpm-1)=AB(IT2-wlc_p%nbpm+1:IT2) ! put end segment type on other end for detail balance +deallocate(Rtemp(3,wlc_p%NB)) +nullify(Utemp(3,wlc_p%NB)) end subroutine