|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 ! vim:fdm=marker: 00002 00003 !------------------------------------------------------------------------------------------- 00010 MODULE modulecreate_V_ext 00011 USE modulegrid 00012 IMPLICIT NONE 00013 00014 00015 CONTAINS 00016 00017 00018 !------------------------------------------------------------------------------------------- 00023 SUBROUTINE write_potentialToFile(NDX,NDY,NDZ,gridx,gridy,gridz,V_ext,filename) 00024 ! {{{ 00025 IMPLICIT NONE 00026 INTEGER,INTENT(IN) :: NDX,NDY,NDZ 00027 REAL*8,INTENT(IN) :: V_ext(NDX*NDY*NDZ) 00028 REAL*8,INTENT(IN) :: gridx(NDX),gridy(NDY),gridz(NDZ) 00029 CHARACTER(150),INTENT(IN) :: filename 00030 INTEGER :: outputunit,ierr 00031 INTEGER :: i,j,k 00032 REAL*8 :: x,y,z 00033 00034 outputunit = 11 00035 ierr = 0 00036 x = 0.d0 00037 y = 0.d0 00038 z = 0.d0 00039 00040 OPEN(UNIT=outputunit,FILE=filename,& 00041 FORM='formatted',ACTION='WRITE',STATUS='REPLACE',IOSTAT=ierr) 00042 00043 openif: IF (ierr == 0) THEN 00044 00045 writeloop: DO k=1,NDZ 00046 DO j=1,NDY 00047 DO i=1,NDX ! format: x,y,z, V(x,y,z) 00048 x=gridx(i) 00049 y=gridy(j) 00050 z=gridz(k) 00051 WRITE(outputunit,*)x,y,z,V_ext(i+NDX*(j-1)+NDX*NDY*(k-1)) 00052 00053 END DO 00054 END DO 00055 END DO writeloop 00056 00057 ELSE openif 00058 WRITE(*,1040) ierr 00059 1040 FORMAT (' ','Error opening file: IOSTAT = ', I6) 00060 END IF openif 00061 00062 CLOSE(outputunit) 00063 00064 00065 ! }}} 00066 END SUBROUTINE write_potentialToFile 00067 00068 00069 !------------------------------------------------------------------------------------------- 00075 SUBROUTINE get_bjjLetterPotential(NDX,NDY,NDZ,gridx,gridy,gridz,Dims,V_ext) 00076 !--- {{{ 00077 IMPLICIT NONE 00078 00079 INTEGER,INTENT(IN) :: NDX,NDY,NDZ,Dims 00080 REAL*8,INTENT(IN) :: gridx(NDX),gridy(NDY),gridz(NDZ) 00081 REAL*8,INTENT(OUT) :: V_ext(NDX*NDY*NDZ) 00082 REAL*8 :: x,y,z 00083 REAL*8 :: x0,xs,a,s0,s2,s3 !for a spline 00084 00085 INTEGER :: i,j,k,m,n 00086 00087 ! this defines a cubic (natural) spline that connects at 00088 ! +-xs to harmonic potentials a(x+-x0)**2. Used in the PRL of Sakmann et al. (2009). 00089 x0 = -2.0d0 00090 xs = -0.5d0 00091 a = 0.5d0 00092 00093 ! for the spline 00094 s0 = 1.d0/3.d0 * (3.d0*a*x0**2 - 2*a*x0*xs) 00095 s2 = (-2.d0*a*x0+a*xs)/xs 00096 s3 = 2.d0/3.d0 * a*x0/xs**2 00097 00098 DO m=1,NDX*NDY*NDZ 00099 00100 CALL get_ijk_from_m(m,NDX,NDY,i,j,k) 00101 x=gridx(i) 00102 y=gridy(j) 00103 z=gridz(k) 00104 00105 IF (x< xs) THEN 00106 V_ext(m)=a*(x-x0)**2 00107 ELSE IF ((x >= xs) .and. (x<=0.d0)) THEN 00108 V_ext(m)=s0 + s2 * x**2 + s3 * x**3 00109 ELSE IF ((x > 0.d0) .and. (x<=-xs)) THEN 00110 V_ext(m)=s0 + s2 * (-x)**2 + s3 * (-x)**3 00111 ELSE IF ((x > -xs)) THEN 00112 V_ext(m)=a*(x+x0)**2 00113 END IF 00114 00115 END DO 00116 00117 !---}}} 00118 END SUBROUTINE get_bjjLetterPotential 00119 00120 !---------------------------------------------------------------------------------------------- 00127 SUBROUTINE get_harmonicPotential(NDX,NDY,NDZ,gridx,gridy,gridz,x0,y0,z0,w_x,w_y,w_z,V_ext) 00128 !--- {{{ 00129 IMPLICIT NONE 00130 00131 INTEGER,INTENT(IN) :: NDX,NDY,NDZ 00132 REAL*8,INTENT(IN) :: gridx(NDX),gridy(NDY),gridz(NDZ),x0,y0,z0,w_x,w_y,w_z 00133 REAL*8,INTENT(OUT) :: V_ext(NDX*NDY*NDZ) 00134 REAL*8 :: x,y,z 00135 00136 INTEGER :: i,j,k,m 00137 00138 DO m=1,NDX*NDY*NDZ 00139 00140 CALL get_ijk_from_m(m,NDX,NDY,i,j,k) 00141 x=gridx(i) 00142 y=gridy(j) 00143 z=gridz(k) 00144 00145 V_ext(m) = 0.5d0 *( w_x**2 * (x-x0)**2 + w_y**2 * (y-y0)**2 + w_z**2 * (z-z0)**2 ) 00146 00147 END DO 00148 00149 !---}}} 00150 END SUBROUTINE get_harmonicPotential 00151 00152 00153 00154 00155 END MODULE modulecreate_V_ext
1.7.3
