hdiff output

r31919/Ackland_metals.f90 2017-02-15 20:30:11.156961897 +0000 r31918/Ackland_metals.f90 2017-02-15 20:30:11.380964909 +0000
  1: !***********************************************************************  1: !***********************************************************************
  2: !SUBROUTINES:                                                          *  2: !SUBROUTINES:                                                          *  
  3: !            BUILD_RHO_SITE                                            *  3: !            BUILD_RHO_SITE                                            *
  4: !            BUILD_V_SITE                                               *  4: !            BUILD_V_SITE                                               *
  5: !***********************************************************************  5: !***********************************************************************
  6: ! FUNCTIONS:  rho_pot,rho_pot_d,rho_pot_dd                             *  6: ! FUNCTIONS:  rho_pot,rho_pot_d,rho_pot_dd                             *
  7: !             Vpot,Vpot_d,Vpot_dd                                       *  7: !             Vpot,Vpot_d,Vpot_dd                                       *             
  8: !             Fembed,Fembed_d,Fembed_dd                                *  8: !             Fembed,Fembed_d,Fembed_dd                                *
  9: !             fcut, fcut_2,fcut_dd                                     *  9: !             fcut, fcut_2,fcut_dd                                     *
 10: !             Mfunc,Mfunc_d, Mfunc_dd,                                 * 10: !             Mfunc,Mfunc_d, Mfunc_dd,                                 *
 11: !             gfunc,gfunc_d,gfunc_dd                                   * 11: !             gfunc,gfunc_d,gfunc_dd                                   *
 12: !             Hfunc                                                    * 12: !             Hfunc                                                    *
 13: !             delta_dirac                                              * 13: !             delta_dirac                                              *
 14: !                                                                      * 14: !                                                                      *             
 15: !****|******************************************************************| 15: !****|******************************************************************|
 16:  16: 
 17: !****|******************************************************************| 17: !****|******************************************************************|
 18: DOUBLE PRECISION FUNCTION rho_pot(ipot,R) 18:         double precision function rho_pot(ipot,R)
 19:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 19:         implicit double precision (a-h,o-z)
 20:   DOUBLE PRECISION Rc 20:         double precision Rc
 21:   INCLUDE 'ackland_sma.h' 21:         include 'ackland_sma.h'  
 22:   INCLUDE 'ackland_mishin_cu.h' 22:         include 'ackland_mishin_cu.h' 
 23:   COMMON /param_cut_off/Rc 23:         COMMON /param_cut_off/Rc        
 24:   INTEGER ipot 24:         integer ipot                
 25:   DOUBLE PRECISION  R 25:         double precision  R
 26:  26:         
 27:   !debug               write(*,*) '================================' 27: !debug               write(*,*) '================================'
 28:   !debug        write(*,*) Rc 28: !debug        write(*,*) Rc
 29:   !debig        stop 29: !debig        stop 
 30:  30:                 
 31:   IF(ipot.EQ.1) THEN 31:         if(ipot.eq.1) then           
 32:      rho_pot=EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta) 32:            rho_pot=exp(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)           
 33:   ELSEIF(ipot.EQ.2) THEN 33:         elseif(ipot.eq.2) then        
 34:      rho_pot=(Ro/R)**(2*q)*fcut(ipot,R,Rc,delta) 34:            rho_pot=(Ro/R)**(2*q)*fcut(ipot,R,Rc,delta)                   
 35:   ELSEIF(ipot.EQ.3) THEN 35:         elseif(ipot.eq.3) then        
 36:      rho_pot=EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta) 36:             rho_pot=exp(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)                   
 37:   ELSEIF(ipot.EQ.4) THEN 37:         elseif(ipot.eq.4) then                        
 38:      rho_pot=gfunc(R,a,beta1,beta2,R03,R04)*fcut(ipot,R,Rc,h) 38:         rho_pot=gfunc(R,a,beta1,beta2,R03,R04)*fcut(ipot,R,Rc,h)        
 39:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6)) THEN 39:         elseif ((ipot.eq.5).or.(ipot.eq.6)) then
 40:      rho_pot = fpsi(R) 40:         rho_pot = fpsi(R)
 41:   ELSEIF ((ipot.EQ.12).OR.(ipot.EQ.13)) THEN 41:         elseif ((ipot.eq.12).or.(ipot.eq.13)) then
 42:      rho_pot =  (   0.77718711248373d0 * (5.6d0-R)**4 & 42:         rho_pot =  (   0.77718711248373d0 * (5.6d0-R)**4 &
 43:           -  0.48102928454986d0 * (5.6d0-R)**5 & 43:                     -  0.48102928454986d0 * (5.6d0-R)**5 &  
 44:           +  0.14501312593993d0 * (5.6d0-R)**6 & 44:                     +  0.14501312593993d0 * (5.6d0-R)**6 &
 45:           -  0.021292226813959d0* (5.6d0-R)**7 & 45:                     -  0.021292226813959d0* (5.6d0-R)**7 &
 46:           +  0.001220921762567d0* (5.6d0-R)**8) * Hfunc(5.6d0-R) 46:                     +  0.001220921762567d0* (5.6d0-R)**8) * Hfunc(5.6d0-R)
 47:   ELSE 47:         else
 48:      WRITE(*,*) 'erreur de ipot' 48:         write(*,*) 'erreur de ipot'        
 49:      rho_pot=0.0D0 49:         
 50:  50:         endif 
 51:   ENDIF 51:            
 52:  52:         Return
 53:   RETURN 53:         End    
 54: END FUNCTION rho_pot 54: 
 55:  55: !****|******************************************************************|        
 56:  56: !****|******************************************************************|        
 57: !****|******************************************************************| 57:        double precision function rho_pot_d(ipot,R)
 58: !****|******************************************************************| 58:         implicit double precision (a-h,o-z)  
 59: DOUBLE PRECISION FUNCTION rho_pot_d(ipot,R) 59:         double precision Rc     
 60:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 60:         include 'ackland_sma.h'  
 61:   DOUBLE PRECISION Rc 61:         include 'ackland_mishin_cu.h' 
 62:   INCLUDE 'ackland_sma.h' 62:         COMMON /param_cut_off/Rc
 63:   INCLUDE 'ackland_mishin_cu.h' 63:         double precision R                
 64:   COMMON /param_cut_off/Rc 64:         integer ipot              
 65:   DOUBLE PRECISION R 65: 
 66:   INTEGER ipot 66: 
 67:  67:         if(ipot.eq.1.) then                                             
 68:  68:       rho_pot_d=-(2*q)/Ro*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)+     &
 69:   IF(ipot.EQ.1.) THEN 69:                EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)
 70:      rho_pot_d=-(2*q)/Ro*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)+     & 70:      
 71:           EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta) 71:      
 72:  72:         elseif(ipot.eq.2) then
 73:  73:     
 74:   ELSEIF(ipot.EQ.2) THEN 74:         rho_pot_d=-2*q/R*(Ro/R)**(2*q)*fcut(ipot,R,Rc,delta)+               &
 75:  75:                (Ro/R)**(2*q)*fcut_d(ipot,R,Rc,delta)        
 76:      rho_pot_d=-2*q/R*(Ro/R)**(2*q)*fcut(ipot,R,Rc,delta)+               & 76:      
 77:           (Ro/R)**(2*q)*fcut_d(ipot,R,Rc,delta) 77:         elseif(ipot.eq.3) then 
 78:  78:                                         
 79:   ELSEIF(ipot.EQ.3) THEN 79:       rho_pot_d=-(2*q)/Ro*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)+     &
 80:  80:                EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)        
 81:      rho_pot_d=-(2*q)/Ro*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)+     & 81:         
 82:           EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta) 82:         elseif(ipot.eq.4) then         
 83:  83:                                
 84:   ELSEIF(ipot.EQ.4) THEN 84:        rho_pot_d=gfunc_d(R,a,beta1,beta2,R03,R04)*fcut(ipot,R,Rc,h)    &
 85:  85:                      +gfunc(R,a,beta1,beta2,R03,R04)*fcut_d(ipot,R,Rc,h)   
 86:      rho_pot_d=gfunc_d(R,a,beta1,beta2,R03,R04)*fcut(ipot,R,Rc,h)    & 86:         elseif ((ipot.eq.5).or.(ipot.eq.6)) then
 87:           +gfunc(R,a,beta1,beta2,R03,R04)*fcut_d(ipot,R,Rc,h) 87:         
 88:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6)) THEN 88:         rho_pot_d=fpsi_d(R)
 89:  89:         elseif ((ipot.eq.12).or.(ipot.eq.13)) then
 90:      rho_pot_d=fpsi_d(R) 90:         rho_pot_d =  (- 0.77718711248373d0 * 4.d0* (5.6d0-R)**3 & 
 91:   ELSEIF ((ipot.EQ.12).OR.(ipot.EQ.13)) THEN 91:                       + 0.48102928454986d0 * 5.d0* (5.6d0-R)**4 &
 92:      rho_pot_d =  (- 0.77718711248373d0 * 4.d0* (5.6d0-R)**3 & 92:                       - 0.14501312593993d0 * 6.d0* (5.6d0-R)**5 &
 93:           + 0.48102928454986d0 * 5.d0* (5.6d0-R)**4 & 93:                       + 0.021292226813959d0* 7.d0* (5.6d0-R)**6 &
 94:           - 0.14501312593993d0 * 6.d0* (5.6d0-R)**5 & 94:                       - 0.001220921762567d0* 8.d0* (5.6d0-R)**7) * Hfunc(5.6d0-R)
 95:           + 0.021292226813959d0* 7.d0* (5.6d0-R)**6 & 95:         else
 96:           - 0.001220921762567d0* 8.d0* (5.6d0-R)**7) * Hfunc(5.6d0-R) 96:         write(*,*) 'erreur de ipot'
 97:   ELSE 97:         endif
 98:      WRITE(*,*) 'erreur de ipot' 98:                 
 99:      rho_pot_d=0.0D0 99:                         
100:   ENDIF100:         Return
101: 101:         End  
102: 102: !****|******************************************************************|        
103:   RETURN103:         double precision  function rho_pot_dd(ipot,R)
104: END FUNCTION rho_pot_d104:         implicit  double precision(a-h,o-z)  
105: !****|******************************************************************|105:         double precision Rc      
106: DOUBLE PRECISION  FUNCTION rho_pot_dd(ipot,R)106:         include 'ackland_sma.h'  
107:   IMPLICIT  DOUBLE PRECISION(a-h,o-z)107:         include 'ackland_mishin_cu.h' 
108:   DOUBLE PRECISION Rc108:         COMMON /param_cut_off/Rc        
109:   INCLUDE 'ackland_sma.h'109:         double precision R                
110:   INCLUDE 'ackland_mishin_cu.h'110:         integer ipot              
111:   COMMON /param_cut_off/Rc111: 
112:   DOUBLE PRECISION R112: 
113:   INTEGER ipot113:        if(ipot.eq.1) then   
114: 114:        rho_pot_dd=(2*q/Ro)**2*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)        &
115: 115:               -2.0d0*(2*q/Ro)*EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)      &
116:   IF(ipot.EQ.1) THEN116:                              +EXP(-2*q*(R/Ro-1.0d0))*fcut_dd(ipot,R,Rc,delta)      
117:      rho_pot_dd=(2*q/Ro)**2*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)        &117:              
118:           -2.0d0*(2*q/Ro)*EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)      &118:         elseif(ipot.eq.2) then
119:           +EXP(-2*q*(R/Ro-1.0d0))*fcut_dd(ipot,R,Rc,delta)119:         rho_pot_dd=2*q*(2*q+1)/R**2*(Ro/R)**(2*q)*fcut(ipot,R,Rc,delta)          &
120: 120:                   -   2*(2*q)/R*(Ro/R)**(2*q)*fcut_d(ipot,R,Rc,delta)            &
121:   ELSEIF(ipot.EQ.2) THEN121:                    + (Ro/R)**(2*q)*fcut_dd(ipot,R,Rc,delta)
122:      rho_pot_dd=2*q*(2*q+1)/R**2*(Ro/R)**(2*q)*fcut(ipot,R,Rc,delta)          &122:        
123:           -   2*(2*q)/R*(Ro/R)**(2*q)*fcut_d(ipot,R,Rc,delta)            &123:         
124:           + (Ro/R)**(2*q)*fcut_dd(ipot,R,Rc,delta)124:         elseif(ipot.eq.3) then
125: 125:        rho_pot_dd=(2*q/Ro)**2*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)         &
126: 126:               -2.0d0*(2*q/Ro)*EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)       &
127:   ELSEIF(ipot.EQ.3) THEN127:                              +EXP(-2*q*(R/Ro-1.0d0))*fcut_dd(ipot,R,Rc,delta)           
128:      rho_pot_dd=(2*q/Ro)**2*EXP(-2*q*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)         &128:                                      
129:           -2.0d0*(2*q/Ro)*EXP(-2*q*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)       &129:         elseif(ipot.eq.4) then         
130:           +EXP(-2*q*(R/Ro-1.0d0))*fcut_dd(ipot,R,Rc,delta)130: !****|******************************************************************|                                       
131: 131:        rho_pot_dd=gfunc_dd(R,a,beta1,beta2,R03,R04)*fcut(ipot,R,Rc,h)     &    
132:   ELSEIF(ipot.EQ.4) THEN132:              +2*gfunc_d(R,a,beta1,beta2,R03,R04)*fcut_d(ipot,R,Rc,h)      &
133:      !****|******************************************************************|133:                 + gfunc(R,a,beta1,beta2,R03,R04)*fcut_dd(ipot,R,Rc,h)  
134:      rho_pot_dd=gfunc_dd(R,a,beta1,beta2,R03,R04)*fcut(ipot,R,Rc,h)     &134: !****|******************************************************************|           
135:           +2*gfunc_d(R,a,beta1,beta2,R03,R04)*fcut_d(ipot,R,Rc,h)      &135:         elseif ((ipot.eq.5).or.(ipot.eq.6)) then
136:           + gfunc(R,a,beta1,beta2,R03,R04)*fcut_dd(ipot,R,Rc,h)136:         rho_pot_dd = fpsi_dd(R)
137:      !****|******************************************************************|137:         
138:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6)) THEN138:         elseif ((ipot.eq.12).or.(ipot.eq.13)) then
139:      rho_pot_dd = fpsi_dd(R)139:         rho_pot_dd =  (  0.77718711248373d0 *12.d0* (5.6d0-R)**2 & 
140: 140:                        - 0.48102928454986d0 *20.d0* (5.6d0-R)**3 &
141:   ELSEIF ((ipot.EQ.12).OR.(ipot.EQ.13)) THEN141:                        + 0.14501312593993d0 *30.d0* (5.6d0-R)**4 &
142:      rho_pot_dd =  (  0.77718711248373d0 *12.d0* (5.6d0-R)**2 &142:                        - 0.021292226813959d0*42.d0* (5.6d0-R)**5 &
143:           - 0.48102928454986d0 *20.d0* (5.6d0-R)**3 &143:                        + 0.001220921762567d0*56.d0* (5.6d0-R)**6) * Hfunc(5.6d0-R)
144:           + 0.14501312593993d0 *30.d0* (5.6d0-R)**4 &144:         else
145:           - 0.021292226813959d0*42.d0* (5.6d0-R)**5 &145:         write(*,*) 'erreur de ipot'
146:           + 0.001220921762567d0*56.d0* (5.6d0-R)**6) * Hfunc(5.6d0-R)146:         endif                
147:   ELSE147:         Return
148:      WRITE(*,*) 'erreur de ipot'148:         End          
149:      rho_pot_dd=0.0D0149:         
150:   ENDIF150: !****|******************************************************************|
151:   RETURN151:         double precision function Vpot(ipot,R)
152: END FUNCTION rho_pot_dd152:         implicit double precision (a-h,o-z)
153: 153:         double precision Rc
154: !****|******************************************************************|154:         include 'ackland_sma.h'  
155: DOUBLE PRECISION FUNCTION Vpot(ipot,R)155:         include 'ackland_mishin_cu.h' 
156:   IMPLICIT DOUBLE PRECISION (a-h,o-z)156:         COMMON /param_cut_off/Rc        
157:   DOUBLE PRECISION Rc157:         integer ipot
158:   INCLUDE 'ackland_sma.h'158:         double precision Mfunc,R        
159:   INCLUDE 'ackland_mishin_cu.h'159:         
160:   COMMON /param_cut_off/Rc160:         if(ipot.eq.1) then           
161:   INTEGER ipot161:            Vpot=exp(-p*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)
162:   DOUBLE PRECISION Mfunc,R162:            Vpot=a_r*Vpot                       
163: 163:         elseif(ipot.eq.2) then        
164:   IF(ipot.EQ.1) THEN164:            Vpot=(Ro/R)**p*fcut(ipot,R,Rc,delta)
165:      Vpot=EXP(-p*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)165:            Vpot=a_r*Vpot                              
166:      Vpot=a_r*Vpot166:         elseif(ipot.eq.3) then                
167:   ELSEIF(ipot.EQ.2) THEN167:            Vpot=(Ro/R)**p*fcut(ipot,R,Rc,delta)        
168:      Vpot=(Ro/R)**p*fcut(ipot,R,Rc,delta)168:            Vpot=a_r*Vpot                                                 
169:      Vpot=a_r*Vpot169:         elseif(ipot.eq.4) then        
170:   ELSEIF(ipot.EQ.3) THEN170:                 
171:      Vpot=(Ro/R)**p*fcut(ipot,R,Rc,delta)171:         Vpot=(E1*Mfunc(R,R01,alpha1)+E2*Mfunc(R,R02,alpha2)+dd)*      &
172:      Vpot=a_r*Vpot172:              fcut(ipot,R,Rc,h)                                        &
173:   ELSEIF(ipot.EQ.4) THEN173:             -S1*Hfunc(Rs1-R)*(Rs1-R)**4                               &
174: 174:             -S2*Hfunc(Rs2-R)*(Rs2-R)**4                               &   
175:      Vpot=(E1*Mfunc(R,R01,alpha1)+E2*Mfunc(R,R02,alpha2)+dd)*      &175:             -S3*Hfunc(Rs3-R)*(Rs3-R)**4        
176:           fcut(ipot,R,Rc,h)                                        &176:         elseif ((ipot.eq.5).or.(ipot.eq.6)) then
177:           -S1*Hfunc(Rs1-R)*(Rs1-R)**4                               &177:         
178:           -S2*Hfunc(Rs2-R)*(Rs2-R)**4                               &178:         Vpot = 0.5d0*fvarphi(R)
179:           -S3*Hfunc(Rs3-R)*(Rs3-R)**4179:         
180:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6)) THEN180:         elseif (ipot.eq.12) then
181: 181:         
182:      Vpot = 0.5d0*fvarphi(R)182:         Vpot = exp(  12.33339230761400d0                 &
183: 183:                    - 10.84732196908600d0*R                 &
184:   ELSEIF (ipot.EQ.12) THEN184:                    +  4.57335244245080d0*R**2            &
185: 185:                    -  0.85266291445935d0*R**3)*            Hfunc(2.3d0-R) * Hfunc0(R-1.d0)  & 
186:      Vpot = EXP(  12.33339230761400d0                 &186:                     
187:           - 10.84732196908600d0*R                 &187:                  +(- 14.261501929757d0*   (3.5d0-R)**4   &
188:           +  4.57335244245080d0*R**2            &188:                    + 15.850036758176d0*   (3.5d0-R)**5   &
189:           -  0.85266291445935d0*R**3)*            Hfunc(2.3d0-R) * Hfunc0(R-1.d0)  &189:                    - 11.325102264291d0*   (3.5d0-R)**6   &
190: 190:                    - 4.0971114831366d0*   (3.5d0-R)**7   &
191:           +(- 14.261501929757d0*   (3.5d0-R)**4   &191:                    + 3.6739378016909d0*   (3.5d0-R)**8  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  &
192:           + 15.850036758176d0*   (3.5d0-R)**5   &192:                  +(  1.3066813393823d0*   (6.0d0-R)**4   &
193:           - 11.325102264291d0*   (3.5d0-R)**6   &193:                    - 0.60542710718094d0*  (6.0d0-R)**5   &
194:           - 4.0971114831366d0*   (3.5d0-R)**7   &194:                    + 1.0055527194350d0 *  (6.0d0-R)**6   &
195:           + 3.6739378016909d0*   (3.5d0-R)**8  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  &195:                    - 0.14918186777562d0*  (6.0d0-R)**7   &
196:           +(  1.3066813393823d0*   (6.0d0-R)**4   &196:                    + 0.032773112059590d0* (6.0d0-R)**8  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)   &
197:           - 0.60542710718094d0*  (6.0d0-R)**5   &197:                  +(  0.011433120304691d0* (7.6d0-R)**4   &
198:           + 1.0055527194350d0 *  (6.0d0-R)**6   &198:                    - 0.021982172508973d0* (7.6d0-R)**5   &
199:           - 0.14918186777562d0*  (6.0d0-R)**7   &199:                    - 0.012542439692607d0* (7.6d0-R)**6   &
200:           + 0.032773112059590d0* (6.0d0-R)**8  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)   &200:                    + 0.025062673874258d0* (7.6d0-R)**7   &
201:           +(  0.011433120304691d0* (7.6d0-R)**4   &201:                    - 0.0075442887837418d0*(7.6d0-R)**8  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0) 
202:           - 0.021982172508973d0* (7.6d0-R)**5   &202: 
203:           - 0.012542439692607d0* (7.6d0-R)**6   &203:         Vpot=0.5d0*Vpot
204:           + 0.025062673874258d0* (7.6d0-R)**7   &204:         elseif (ipot.eq.13) then
205:           - 0.0075442887837418d0*(7.6d0-R)**8  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0)205:         
206: 206:         Vpot = exp(  12.8822300381920d0                         &
207:      Vpot=0.5d0*Vpot207:                    - 12.1838501578140d0*R                       &
208:   ELSEIF (ipot.EQ.13) THEN208:                    +  5.5998956281737d0*R**2                     &
209: 209:                    -  1.0915156420318d0*R**3)               *  Hfunc(2.3d0-R)*Hfunc0(R-1.d0)  & 
210:      Vpot = EXP(  12.8822300381920d0                         &210:                    
211:           - 12.1838501578140d0*R                       &211:                  +(   8.4670497139946d0*     (3.5d0-R)**4    &
212:           +  5.5998956281737d0*R**2                     &212:                    - 46.183472786003d0*      (3.5d0-R)**5    &
213:           -  1.0915156420318d0*R**3)               *  Hfunc(2.3d0-R)*Hfunc0(R-1.d0)  &213:                    + 79.633499844770d0*      (3.5d0-R)**6    &
214: 214:                    - 64.847634731465d0*      (3.5d0-R)**7    &
215:           +(   8.4670497139946d0*     (3.5d0-R)**4    &215:                    + 19.454623850774d0*      (3.5d0-R)**8  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  & 
216:           - 46.183472786003d0*      (3.5d0-R)**5    &216:                    
217:           + 79.633499844770d0*      (3.5d0-R)**6    &217:                  +(-  0.097845860135187d0*   (6.0d0-R)**4    &
218:           - 64.847634731465d0*      (3.5d0-R)**7    &218:                    -  0.47537134413743d0*    (6.0d0-R)**5    &
219:           + 19.454623850774d0*      (3.5d0-R)**8  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  &219:                    -  0.00096806164225329d0* (6.0d0-R)**6    &
220: 220:                    -  0.16355187497617d0*    (6.0d0-R)**7    &
221:           +(-  0.097845860135187d0*   (6.0d0-R)**4    &221:                    -  0.00090914903435333d0* (6.0d0-R)**8  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)   &   
222:           -  0.47537134413743d0*    (6.0d0-R)**5    &222:                      
223:           -  0.00096806164225329d0* (6.0d0-R)**6    &223:                  +(-  0.022038480751134d0*   (7.6d0-R)**4    &
224:           -  0.16355187497617d0*    (6.0d0-R)**7    &224:                    -  0.060955465943384d0*   (7.6d0-R)**5    &
225:           -  0.00090914903435333d0* (6.0d0-R)**8  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)   &225:                    +  0.11573689045653d0*    (7.6d0-R)**6    &
226: 226:                    -  0.062697675088029d0*   (7.6d0-R)**7    &
227:           +(-  0.022038480751134d0*   (7.6d0-R)**4    &227:                    +  0.011273545085049d0*   (7.6d0-R)**8  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0) 
228:           -  0.060955465943384d0*   (7.6d0-R)**5    &228: 
229:           +  0.11573689045653d0*    (7.6d0-R)**6    &229:         Vpot=0.5d0*Vpot
230:           -  0.062697675088029d0*   (7.6d0-R)**7    &230:         else
231:           +  0.011273545085049d0*   (7.6d0-R)**8  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0)231:         write(*,*) 'erreur de ipot'        
232: 232:         
233:      Vpot=0.5d0*Vpot233:         endif    
234:   ELSE234:            
235:      WRITE(*,*) 'erreur de ipot'235:         Return
236:      Vpot=0.0D0236:         End    
237: 237: !****|******************************************************************|
238:   ENDIF238:         double precision function Vpot_d(ipot,R)
239: 239:         implicit double precision (a-h,o-z)        
240:   RETURN240:         include 'ackland_sma.h'  
241: END FUNCTION Vpot241:         include 'ackland_mishin_cu.h'        
242: !****|******************************************************************|242:         double precision Rc
243: DOUBLE PRECISION FUNCTION Vpot_d(ipot,R)243:         COMMON /param_cut_off/Rc         
244:   IMPLICIT DOUBLE PRECISION (a-h,o-z)244:         double precision Mfunc,Mfunc_d,R                
245:   INCLUDE 'ackland_sma.h'245:         integer ipot              
246:   INCLUDE 'ackland_mishin_cu.h'246: 
247:   DOUBLE PRECISION Rc247: 
248:   COMMON /param_cut_off/Rc248:         if(ipot.eq.1) then
249:   DOUBLE PRECISION Mfunc,Mfunc_d,R249:                                               
250:   INTEGER ipot250:         Vpot_d=-p/Ro*DEXP(-p*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)+    &
251: 251:                DEXP(-p*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)
252: 252:         Vpot_d=a_r*Vpot_d     
253:   IF(ipot.EQ.1) THEN253:      
254: 254:         elseif(ipot.eq.2) then
255:      Vpot_d=-p/Ro*DEXP(-p*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)+    &255:     
256:           DEXP(-p*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)256:         Vpot_d=-p/R*(Ro/R)**p*fcut(ipot,R,Rc,delta)+                 &
257:      Vpot_d=a_r*Vpot_d257:                (Ro/R)**p*fcut_d(ipot,R,Rc,delta)
258: 258:         Vpot_d=a_r*Vpot_d     
259:   ELSEIF(ipot.EQ.2) THEN259:         elseif(ipot.eq.3) then
260: 260:     
261:      Vpot_d=-p/R*(Ro/R)**p*fcut(ipot,R,Rc,delta)+                 &261:         Vpot_d=-p/R*(Ro/R)**p*fcut(ipot,R,Rc,delta)+                 &
262:           (Ro/R)**p*fcut_d(ipot,R,Rc,delta)262:                (Ro/R)**p*fcut_d(ipot,R,Rc,delta)                
263:      Vpot_d=a_r*Vpot_d263:         Vpot_d=a_r*Vpot_d    
264:   ELSEIF(ipot.EQ.3) THEN264:         elseif(ipot.eq.4) then 
265: 265: 
266:      Vpot_d=-p/R*(Ro/R)**p*fcut(ipot,R,Rc,delta)+                 &266:       Vpot_d=(E1*Mfunc_d(R,R01,alpha1)+E2*Mfunc_d(R,R02,alpha2))*    &
267:           (Ro/R)**p*fcut_d(ipot,R,Rc,delta)267:             fcut(ipot,R,Rc,h)                                        &
268:      Vpot_d=a_r*Vpot_d268:             +4*S1*Hfunc(Rs1-R)*(Rs1-R)**3                            &
269:   ELSEIF(ipot.EQ.4) THEN269:             +4*S2*Hfunc(Rs2-R)*(Rs2-R)**3                            &    
270: 270:             +4*S3*Hfunc(Rs3-R)*(Rs3-R)**3                            &
271:      Vpot_d=(E1*Mfunc_d(R,R01,alpha1)+E2*Mfunc_d(R,R02,alpha2))*    &271:       +    (E1*Mfunc(R,R01,alpha1)+E2*Mfunc(R,R02,alpha2)+dd)*       &
272:           fcut(ipot,R,Rc,h)                                        &272:             fcut_d(ipot,R,Rc,h)                                 
273:           +4*S1*Hfunc(Rs1-R)*(Rs1-R)**3                            &273:         elseif ((ipot.eq.5).or.(ipot.eq.6)) then
274:           +4*S2*Hfunc(Rs2-R)*(Rs2-R)**3                            &274:         
275:           +4*S3*Hfunc(Rs3-R)*(Rs3-R)**3                            &275:         Vpot_d = 0.5d0*fvarphi_d(R)
276:           +    (E1*Mfunc(R,R01,alpha1)+E2*Mfunc(R,R02,alpha2)+dd)*       &276:      
277:           fcut_d(ipot,R,Rc,h)277:         elseif (ipot.eq.12) then
278:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6)) THEN278:         
279: 279:         Vpot_d= exp( 12.333392307614d0           -  10.847321969086d0*R            &
280:      Vpot_d = 0.5d0*fvarphi_d(R)280:                     +4.5733524424508d0*R**2 -  0.85266291445935d0*R**3)     &
281: 281:                 *( -  10.847321969086d0           +   4.5733524424508d0*2.d0*R     &
282:   ELSEIF (ipot.EQ.12) THEN282:                    -  0.85266291445935d0*3.d0*R**2)                         &   
283: 283:                     * Hfunc(2.3d0-R) * Hfunc(R-1.d0)                        &  
284:      Vpot_d= EXP( 12.333392307614d0           -  10.847321969086d0*R            &284:                  +(  14.261501929757d0*   4.d0*  (3.5d0-R)**3   &
285:           +4.5733524424508d0*R**2 -  0.85266291445935d0*R**3)     &285:                    - 15.850036758176d0*   5.d0*  (3.5d0-R)**4   &
286:           *( -  10.847321969086d0           +   4.5733524424508d0*2.d0*R     &286:                    + 11.325102264291d0*   6.d0*  (3.5d0-R)**5   &
287:           -  0.85266291445935d0*3.d0*R**2)                         &287:                    + 4.0971114831366d0*   7.d0*  (3.5d0-R)**6   &
288:           * Hfunc(2.3d0-R) * Hfunc(R-1.d0)                        &288:                    - 3.6739378016909d0*   8.d0*  (3.5d0-R)**7  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  &
289:           +(  14.261501929757d0*   4.d0*  (3.5d0-R)**3   &289:                  +(- 1.3066813393823d0*   4.d0*  (6.0d0-R)**3   &
290:           - 15.850036758176d0*   5.d0*  (3.5d0-R)**4   &290:                    + 0.60542710718094d0*  5.d0*  (6.0d0-R)**4   &
291:           + 11.325102264291d0*   6.d0*  (3.5d0-R)**5   &291:                    - 1.0055527194350d0*   6.d0*  (6.0d0-R)**5   &
292:           + 4.0971114831366d0*   7.d0*  (3.5d0-R)**6   &292:                    + 0.14918186777562d0*  7.d0*  (6.0d0-R)**6   &
293:           - 3.6739378016909d0*   8.d0*  (3.5d0-R)**7  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  &293:                    - 0.032773112059590*   8.d0*  (6.0d0-R)**7  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)        &
294:           +(- 1.3066813393823d0*   4.d0*  (6.0d0-R)**3   &294:                  +(- 0.011433120304691d0* 4.d0*  (7.6d0-R)**3   &
295:           + 0.60542710718094d0*  5.d0*  (6.0d0-R)**4   &295:                    + 0.021982172508973d0* 5.d0*  (7.6d0-R)**4   &
296:           - 1.0055527194350d0*   6.d0*  (6.0d0-R)**5   &296:                    + 0.012542439692607d0* 6.d0*  (7.6d0-R)**5   &
297:           + 0.14918186777562d0*  7.d0*  (6.0d0-R)**6   &297:                    - 0.025062673874258d0* 7.d0*  (7.6d0-R)**6   &
298:           - 0.032773112059590*   8.d0*  (6.0d0-R)**7  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)        &298:                    + 0.0075442887837418d0*8.d0*  (7.6d0-R)**7  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0) 
299:           +(- 0.011433120304691d0* 4.d0*  (7.6d0-R)**3   &299: 
300:           + 0.021982172508973d0* 5.d0*  (7.6d0-R)**4   &300:         Vpot_d=0.5d0*Vpot_d
301:           + 0.012542439692607d0* 6.d0*  (7.6d0-R)**5   &301:         elseif (ipot.eq.13) then
302:           - 0.025062673874258d0* 7.d0*  (7.6d0-R)**6   &302:         
303:           + 0.0075442887837418d0*8.d0*  (7.6d0-R)**7  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0)303:         Vpot_d = exp( 12.882230038192d0      -   12.183850157814*R               &
304: 304:                    +5.5998956281737d0*R**2 -    1.0915156420318d0*R**3)        &
305:      Vpot_d=0.5d0*Vpot_d305:                  *(-12.183850157814        +    5.5998956281737d0*2.d0*R       &
306:   ELSEIF (ipot.EQ.13) THEN306:                    -1.0915156420318d0*3.d0*R**2)                               &
307: 307:                    * Hfunc(2.3d0-R) * Hfunc0(R-1.d0)                            &
308:      Vpot_d = EXP( 12.882230038192d0      -   12.183850157814*R               &308:                  +(- 8.4670497139946d0*     4.d0*  (3.5d0-R)**3    &
309:           +5.5998956281737d0*R**2 -    1.0915156420318d0*R**3)        &309:                    + 46.183472786003d0*     5.d0*  (3.5d0-R)**4    &
310:           *(-12.183850157814        +    5.5998956281737d0*2.d0*R       &310:                    - 79.633499844770d0*     6.d0*  (3.5d0-R)**5    &
311:           -1.0915156420318d0*3.d0*R**2)                               &311:                    + 64.847634731465d0*     7.d0*  (3.5d0-R)**6    &
312:           * Hfunc(2.3d0-R) * Hfunc0(R-1.d0)                            &312:                    - 19.454623850774d0*     8.d0*  (3.5d0-R)**7  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  & 
313:           +(- 8.4670497139946d0*     4.d0*  (3.5d0-R)**3    &313:                  +(+ 0.097845860135187d0*   4.d0*  (6.0d0-R)**3    &
314:           + 46.183472786003d0*     5.d0*  (3.5d0-R)**4    &314:                    + 0.47537134413743d0*    5.d0*  (6.0d0-R)**4    &
315:           - 79.633499844770d0*     6.d0*  (3.5d0-R)**5    &315:                    + 0.00096806164225329d0* 6.d0*  (6.0d0-R)**5    &
316:           + 64.847634731465d0*     7.d0*  (3.5d0-R)**6    &316:                    + 0.16355187497617d0*    7.d0*  (6.0d0-R)**6    &
317:           - 19.454623850774d0*     8.d0*  (3.5d0-R)**7  ) * Hfunc(3.5d0-R)*Hfunc(R-2.3d0)  &317:                    + 0.00090914903435333d0* 8.d0*  (6.0d0-R)**7  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)   &          
318:           +(+ 0.097845860135187d0*   4.d0*  (6.0d0-R)**3    &318:                  +(+ 0.022038480751134d0*   4.d0*  (7.6d0-R)**3    &
319:           + 0.47537134413743d0*    5.d0*  (6.0d0-R)**4    &319:                    + 0.060955465943384d0*   5.d0*  (7.6d0-R)**4    &
320:           + 0.00096806164225329d0* 6.d0*  (6.0d0-R)**5    &320:                    - 0.11573689045653d0*    6.d0*  (7.6d0-R)**5    &
321:           + 0.16355187497617d0*    7.d0*  (6.0d0-R)**6    &321:                    + 0.062697675088029d0*   7.d0*  (7.6d0-R)**6    &
322:           + 0.00090914903435333d0* 8.d0*  (6.0d0-R)**7  ) * Hfunc(6.d0-R)*Hfunc(R-2.3d0)   &322:                    - 0.011273545085049d0*   8.d0*  (7.6d0-R)**7  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0) 
323:           +(+ 0.022038480751134d0*   4.d0*  (7.6d0-R)**3    &323:         Vpot_d=0.5d0*Vpot_d
324:           + 0.060955465943384d0*   5.d0*  (7.6d0-R)**4    &324: 
325:           - 0.11573689045653d0*    6.d0*  (7.6d0-R)**5    &325:         else
326:           + 0.062697675088029d0*   7.d0*  (7.6d0-R)**6    &326:         write(*,*) 'erreur de ipot'
327:           - 0.011273545085049d0*   8.d0*  (7.6d0-R)**7  ) * Hfunc(7.6d0-R)*Hfunc(R-2.3d0)327:         endif                                 
328:      Vpot_d=0.5d0*Vpot_d328:     
329: 329:         Return
330:   ELSE330:         End    
331:      WRITE(*,*) 'erreur de ipot'331:         
332:      Vpot_d=0.0D0332: !****|******************************************************************|
333:   ENDIF333:         double precision function Vpot_dd(ipot,R)
334: 334:         implicit double precision (a-h,o-z)        
335:   RETURN335:         include 'ackland_sma.h'  
336: END FUNCTION Vpot_d336:         include 'ackland_mishin_cu.h' 
337: 337:         double precision Rc
338: !****|******************************************************************|338:         COMMON /param_cut_off/Rc        
339: DOUBLE PRECISION FUNCTION Vpot_dd(ipot,R)339:         double precision Mfunc,Mfunc_d,Mfunc_dd,R                
340:   IMPLICIT DOUBLE PRECISION (a-h,o-z)340:         integer ipot              
341:   INCLUDE 'ackland_sma.h' 
342:   INCLUDE 'ackland_mishin_cu.h' 
343:   DOUBLE PRECISION Rc 
344:   COMMON /param_cut_off/Rc 
345:   DOUBLE PRECISION Mfunc,Mfunc_d,Mfunc_dd,R 
346:   INTEGER ipot 
347: 341: 
348: 342: 
349:   IF(ipot.EQ.1) THEN343:         if(ipot.eq.1) then   
350:      Vpot_dd=(p/Ro)**2*EXP(-p*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)      &344:        Vpot_dd=(p/Ro)**2*EXP(-p*(R/Ro-1.0d0))*fcut(ipot,R,Rc,delta)      &
351:           -2.0d0*(p/Ro)*EXP(-p*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)     &345:           -2.0d0*(p/Ro)*EXP(-p*(R/Ro-1.0d0))*fcut_d(ipot,R,Rc,delta)     &
352:           +EXP(-p*(R/Ro-1.0d0))*fcut_dd(ipot,R,Rc,delta)346:               +EXP(-p*(R/Ro-1.0d0))*fcut_dd(ipot,R,Rc,delta)
353:      Vpot_dd=a_r*Vpot_dd347:          Vpot_dd=a_r*Vpot_dd   
 348:              
 349:         elseif(ipot.eq.2) then
 350:         Vpot_dd=p*(p+1)/R**2*(Ro/R)**p*fcut(ipot,R,Rc,delta)          &
 351:                       -2*p/R*(Ro/R)**p*fcut_d(ipot,R,Rc,delta)        &
 352:                             +(Ro/R)**p*fcut_dd(ipot,R,Rc,delta)       
 353:         Vpot_dd=a_r*Vpot_dd  
 354:         
 355:         elseif(ipot.eq.3) then
 356:        Vpot_dd=p*(p+1)/R**2*(Ro/R)**p*fcut(ipot,R,Rc,delta)           &
 357:                       -2*p/R*(Ro/R)**p*fcut_d(ipot,R,Rc,delta)        &
 358:                             +(Ro/R)**p*fcut_dd(ipot,R,Rc,delta)       
 359:          Vpot_dd=a_r*Vpot_dd  
 360:                             
 361:         elseif(ipot.eq.4) then 
 362: 
 363: !****|******************************************************************|                                       
 364:       Vpot_dd=(E1*Mfunc_dd(R,R01,alpha1)+E2*Mfunc_dd(R,R02,alpha2))*fcut(ipot,R,Rc,h)       &
 365:        +    2*(E1*Mfunc_d(R,R01,alpha1)+E2*Mfunc_d(R,R02,alpha2))*fcut_d(ipot,R,Rc,h)       &
 366:        +     (E1*Mfunc(R,R01,alpha1)+E2*Mfunc(R,R02,alpha2)+dd)*fcut_dd(ipot,R,Rc,h)        &   
 367:        -      12*S1*Hfunc(Rs1-R)*(Rs1-R)**2                                                        &
 368:        -      12*S2*Hfunc(Rs2-R)*(Rs2-R)**2                                                   &
 369:        -      12*S3*Hfunc(Rs3-R)*(Rs3-R)**2                                                        
 370:                       
 371: !****|******************************************************************|           
 372:          elseif ((ipot.eq.5).or.(ipot.eq.6)) then
 373:         
 374:         Vpot_dd = 0.5d0*fvarphi_dd(R)
 375:         
 376:         else
 377:         write(*,*) 'erreur de ipot'
 378:         endif                                 
 379:     
 380:         Return
 381:         End                    
 382:         
 383:           
 384: !****|******************************************************************|        
 385:         double precision  function fcut(ipot,R,Rc,width)
 386:         implicit double precision (a-h,o-z)        
 387:         double precision R,Rc,width,test,x
 388:         
 389:         if(ipot.eq.1.or.ipot.eq.2.or.ipot.eq.3) then
 390:            test=(R-Rc)/width
 391:            if(test.lt.-100.d0)  then
 392:            fcut=1.0d0        
 393:            elseif(test.gt.100.d0) then
 394:            fcut=0.0d0        
 395:            else        
 396:            x=test        
 397:            fcut=1.0d0/(1.0d0+EXP(x))
 398:            endif
 399:         elseif(ipot.eq.4) then           
 400:             test=(R-Rc)
 401:            if(test.gt.0.d0)  then
 402:            fcut=0.0d0
 403:            elseif(test.le.0.d0) then
 404:            x=test/width           
 405:            fcut= x**4/(1+x**4)                      
 406:            endif
 407:         elseif ((ipot.eq.5).or.(ipot.eq.6).or.(ipot.eq.12).or.(ipot.eq.13)) then
 408:         fcut=1.d0                      
 409:         endif
 410:                                    
 411:         Return
 412:         End           
 413:         
 414: !****|******************************************************************|        
 415:         double precision  function fcut_d(ipot,R,Rc,width)
 416:         implicit double precision (a-h,o-z)        
 417:         double precision R,Rc,width,test,x
 418:         
 419:         if(ipot.eq.1.or.ipot.eq.2.or.ipot.eq.3) then        
 420:           test=(R-Rc)/width
 421:            if(test.lt.-100.d0)  then
 422:              fcut_d=0.0d0        
 423:            elseif(test.gt.100.d0) then
 424:              fcut_d=0.0d0                
 425:             else
 426:              x=test                
 427:              fcut_d=-EXP(x)/(1+EXP(x))**2
 428:              fcut_d=fcut_d/width          
 429:            endif                
 430:         elseif(ipot.eq.4) then        
 431:             test=(R-Rc)
 432:            if(test.gt.0.d0)  then
 433:            fcut_d=0.0d0        
 434:            elseif(test.le.0.d0) then
 435:            x=test/width
 436:            fcut_d= 4*x**3/(1+x**4)**2           
 437:            fcut_d=fcut_d/width               
 438:            endif                
 439:         elseif ((ipot.eq.5).or.(ipot.eq.6).or.(ipot.eq.12).or.(ipot.eq.13)) then
 440:         fcut_d=1.d0                      
 441:         endif        
 442:                 
 443:         Return
 444:         End           
 445:                            
 446: !****|******************************************************************|        
 447:         double precision function fcut_dd(ipot,R,Rc,width)
 448:         implicit  double precision(a-h,o-z)        
 449:         double precision R,Rc,width,test,x
 450:         
 451:         if(ipot.eq.1.or.ipot.eq.2.or.ipot.eq.3) then        
 452:           test=(R-Rc)/width
 453:            if(test.lt.-100.d0)  then
 454:              fcut_dd=0.0d0        
 455:            elseif(test.gt.100.d0) then
 456:              fcut_dd=0.0d0                
 457:             else
 458:              x=test                
 459:              fcut_dd= -EXP(x)/( (1.0d0+EXP(x) ) )**2 +2.0d0*EXP(2*x)/(1.0d0+EXP(x) )**3    
 460:            fcut_dd=fcut_dd/width**2     
 461:            endif                
 462:         elseif(ipot.eq.4) then        
 463:             test=(R-Rc)
 464:            if(test.gt.0.d0)  then
 465:            fcut_dd=0.0d0        
 466:            elseif(test.le.0.d0) then
 467:            x=test/width           
 468:            fcut_dd=12*x**2/(1+x**4)**2-32*x**6/(1+x**4)**3
 469:            fcut_dd=fcut_dd/width**2                          
 470:            endif                
 471:         elseif ((ipot.eq.5).or.(ipot.eq.6).or.(ipot.eq.12).or.(ipot.eq.13)) then
 472:         fcut_dd=1.d0                      
 473:         
 474:         endif        
 475:                 
 476:         Return
 477:         End          
 478:                                              
 479: !****|******************************************************************|
 480:         double precision function Fembed(ipot,x)
 481:         USE COMMONS, ONLY : ACK1, ACK2
 482:         implicit double precision (a-h,o-z)        
 483:         double precision x
 484:         integer ipot
 485:         double precision Rc
 486:         COMMON /param_cut_off/Rc        
 487:         include 'ackland_sma.h'  
 488:         include 'ackland_mishin_cu.h'
 489:         include 'ackland_mendelev_fe.h'        
 490: 
 491:                 
 492:           if(ipot.eq.1.or.ipot.eq.2.or.ipot.eq.3) then        
 493:              Fembed=x**(alpha)
 494:              Fembed=-b_a*Fembed          
 495:                 
 496:           elseif(ipot.eq.4) then             
 497:           
 498:           if(x.le.1) then                    
 499:            Fembed=F0+0.5*F2*(x-1)**2+q1*(x-1)**3+q2*(x-1)**4 +q3*(x-1)**5 +q4*(x-1)**6 
 500:         
 501:           elseif(x.gt.1) then                            
 502:            Fembed=(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)/(1+qq2*(x-1)**3)               
 503:           endif                     
 504:                     
 505:           elseif(ipot.eq.5) then
 506:           
 507:           Fembed = -dsqrt(x) + aphi*x*x + aphi2*x*x*x*x  
 508:           
 509:           elseif(ipot.eq.6) then
 510:            if (x.lt.1.d0) then
 511:              Fembed = -aphi*dsqrt(x) - aphi2*(1-dsqrt(x))*dlog(2.d0-x)/dlog(2.d0) 
 512:              else
 513:              Fembed = -aphi*dsqrt(x) 
 514:            endif
 515: 
 516:           elseif(ipot.eq.12) then
 517:            Fembed = -     dsqrt(x)                                         &
 518:                     - 1.9162462126235d0*1.d-7*(x-60.d0)**4*Hfunc(x-60.d0) &
 519:                     + 4.6418727035037d0*1.d-7*(x-70.d0)**4*Hfunc(x-70.d0) &
 520:                     + 6.6448294272955d0*1.d-7*(x-80.d0)**4*Hfunc(x-80.d0) &
 521:                     - 2.0680252960229d0*1.d-6*(x-85.d0)**4*Hfunc(x-85.d0) &
 522:                     + 1.1387131464983d0*1.d-6*(x-90.d0)**4*Hfunc(x-90.d0)
 523:           elseif(ipot.eq.13) then
 524:            Fembed = - dsqrt(x)                                         &
 525:                     + 3.2283012597866d0*1.d-7*(x-60.d0)**4*Hfunc(x-60.d0) &                
 526:                     - 1.1552813894483d0*1.d-6*(x-70.d0)**4*Hfunc(x-70.d0) &
 527:                     + 2.3747280268355d0*1.d-6*(x-80.d0)**4*Hfunc(x-80.d0) &
 528:                     - 2.0379550826523d0*1.d-6*(x-85.d0)**4*Hfunc(x-85.d0) & 
 529:                     + 4.9758343293936d0*1.d-7*(x-90.d0)**4*Hfunc(x-90.d0)
 530: 
 531: 
 532:          end if 
 533:         return
 534:         end
 535: !****|******************************************************************|
 536:         double precision function Fembed_d(ipot,x)
 537:         USE COMMONS, ONLY : ACK1, ACK2
 538:         implicit double precision (a-h,o-z)        
 539:         double precision x
 540:         double precision Rc
 541:         COMMON /param_cut_off/Rc        
 542:         include 'ackland_sma.h'  
 543:         include 'ackland_mishin_cu.h'
 544:         include 'ackland_mendelev_fe.h'
 545:         
 546:           
 547:           if(ipot.eq.1.or.ipot.eq.2.or.ipot.eq.3) then        
 548:              Fembed_d=alpha*x**(alpha-1.0d0)
 549:              Fembed_d=-b_a*Fembed_d                       
 550:           elseif(ipot.eq.4) then                       
 551:           if(x.le.1) then                    
 552:            Fembed_d=F2*(x-1)+3*q1*(x-1)**2+4*q2*(x-1)**3+5*q3*(x-1)**4+6*q4*(x-1)**5 
 553:         
 554:           elseif(x.gt.1) then                            
 555:          Fembed_d=(F2*(x-1)+3*q1*(x-1)**2+4*qq1*(x-1)**3)/(1+qq2*(x-1)**3)       &
 556:                   -3*qq2*(x-1)**2*(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)  &
 557:                   /(1+qq2*(x-1)**3)**2      
 558:           
 559:                          
 560:            endif                               
 561: 
 562:           elseif (ipot.eq.5) then
 563:           
 564:           Fembed_d = -0.5d0/dsqrt(x) + 2.d0*aphi*x  &
 565:                       + 4.d0*aphi2*x*x*x
 566: 
 567: 
 568:           elseif(ipot.eq.6) then
 569:           if (x.lt.1.d0) then
 570:             Fembed_d = -aphi*0.5d0/dsqrt(x)                   &
 571:                        + aphi2*(0.5d0*dlog(2.d0-x)/dsqrt(x)  &
 572:                        + (1-dsqrt(x))/(2.d0-x))/dlog(2.d0) 
 573:            else
 574:             Fembed_d = -0.5d0*aphi/dsqrt(x)     
 575:           end if             
 576:             
 577:           elseif(ipot.eq.12) then
 578:            Fembed_d = - 0.5d0 / dsqrt(x)                                 &
 579:                     + 4.d0*(                                           &
 580:                     - 1.9162462126235d0*1.d-7*(x-60.d0)**3*Hfunc(x-60.d0) &
 581:                     + 4.6418727035037d0*1.d-7*(x-70.d0)**3*Hfunc(x-70.d0) &
 582:                     + 6.6448294272955d0*1.d-7*(x-80.d0)**3*Hfunc(x-80.d0) &
 583:                     - 2.0680252960229d0*1.d-6*(x-85.d0)**3*Hfunc(x-85.d0) &
 584:                     + 1.1387131464983d0*1.d-6*(x-90.d0)**3*Hfunc(x-90.d0) &
 585:                            )
 586:           elseif(ipot.eq.13) then
 587:            Fembed_d = - 0.5d0 / dsqrt(x)                                 &
 588:                     + 4.d0*(                                           &
 589:                     + 3.2283012597866d0*1.d-7*(x-60.d0)**3*Hfunc(x-60.d0) &                
 590:                     - 1.1552813894483d0*1.d-6*(x-70.d0)**3*Hfunc(x-70.d0) &
 591:                     + 2.3747280268355d0*1.d-6*(x-80.d0)**3*Hfunc(x-80.d0) &
 592:                     - 2.0379550826523d0*1.d-6*(x-85.d0)**3*Hfunc(x-85.d0) & 
 593:                     + 4.9758343293936d0*1.d-7*(x-90.d0)**3*Hfunc(x-90.d0) &
 594:                            )
 595:         
 596:         
 597:         endif
 598: 
 599:         return
 600:         end
 601:         
 602: !****|******************************************************************|        
 603:         double precision function Fembed_dd(ipot,x)
 604:         USE COMMONS, ONLY : ACK1, ACK2
 605:         implicit double precision (a-h,o-z)        
 606:         include 'ackland_sma.h'  
 607:         include 'ackland_mishin_cu.h'
 608:         double precision Rc
 609:         COMMON /param_cut_off/Rc        
 610:         double precision x,denom
 611:         include 'ackland_mendelev_fe.h'
 612:           
 613:           if(ipot.eq.1.or.ipot.eq.2.or.ipot.eq.3) then        
 614:              Fembed_dd=alpha*(alpha-1)*x**(alpha-2.0d0)
 615:              Fembed_dd=-b_a*Fembed_dd                       
 616:           elseif(ipot.eq.4) then                       
 617:           if(x<1) then
 618:            Fembed_dd=F2+6*q1*(x-1)+12*q2*(x-1)**2+20*q3*(x-1)**3+30*q4*(x-1)**4 
 619:         
 620:           elseif(x>1) then
 621:          denom=1+qq2*(x-1)**3
 622:                                                
 623:          Fembed_dd=(F2+6*q1*(x-1)+12*qq1*(x-1)**2)/denom                              &
 624:          -6*qq2*(x-1)**2*(F2*(x-1)+3*q1*(x-1)**2+4*qq1*(x-1)**3)/denom**2             &
 625:          -6*qq2*(x-1)*(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)/denom**2          &
 626:           +18*qq2**2*(x-1)**4*(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)/denom**3  
 627:                                    
 628:            endif                               
 629:                     
 630:           elseif (ipot.eq.5) then
 631:             Fembed_dd = 0.25d0/dsqrt(x**3) + 2.d0*aphi  &
 632:                       + 12.d0*aphi2*x*x
 633:                       
 634:           elseif (ipot.eq.6) then
 635:            if (x.lt.1.d0) then
 636:              Fembed_dd = 0.25d0*aphi/dsqrt(x**3)  + aphi2/dlog(2.d0)*   &
 637:                      (                                              &
 638:                        -0.25d0*dlog(2.d0-x)/dsqrt(x**3)             &
 639:                        -1.0/((2.d0-x)*dsqrt(x))                     &
 640:                        +(1.0-dsqrt(x))/(2.d0-x)**2                  &
 641:                       ) 
 642:              else
 643:               Fembed_dd = 0.25d0*aphi/dsqrt(x**3)     
 644:             end if 
 645:          
 646:         endif
 647: 
 648:         return
 649:         end
 650: !****|******************************************************************|
 651:         double precision function delta_dirac(i,j)
 652:         integer i,j
 653: 
 654:          delta_dirac=0.0d0
 655:         if(i.eq.j) then
 656:          delta_dirac=1.0d0
 657:         endif
 658:         return
 659:         end
 660: !****|******************************************************************|        
 661:         double precision function Hfunc(x)
 662:         double precision x
 663:         
 664:         if(x.lt.0.0) then
 665:           Hfunc=0.0d0
 666:         elseif(x.ge.0.0) then          
 667:           Hfunc=1.0d0         
 668:          endif
 669:         return
 670:         end         
 671: !****|******************************************************************|        
 672:         double precision function Hfunc0(x)
 673:         double precision x
 674:         
 675:         if(x.le.0.0) then
 676:           Hfunc0=0.0d0
 677:         elseif(x.gt.0.0) then          
 678:           Hfunc0=1.0d0         
 679:          endif
 680:         return
 681:         end         
 682: !****|******************************************************************|
 683:         double precision function Mfunc(R,R0,alpha)
 684:         double precision R,R0,alpha
 685:         
 686:         Mfunc=exp(-2*alpha*(R-R0))-2*exp(-alpha*(R-R0))
 687:                 
 688:         return
 689:         end         
 690: !****|******************************************************************|
 691:         double precision function Mfunc_d(R,R0,alpha)
 692:         double precision R,R0,alpha
 693:         
 694:        Mfunc_d=-2*alpha*exp(-2*alpha*(R-R0))+2*alpha*exp(-alpha*(R-R0))
 695:                 
 696:         return
 697:         end
 698: !****|******************************************************************|        
 699:         double precision function Mfunc_dd(R,R0,alpha)
 700:         double precision R,R0,alpha
 701:         
 702:          Mfunc_dd=+4*alpha**2*exp(-2*alpha*(R-R0))-2*alpha**2*exp(-alpha*(R-R0))
 703:                 
 704:         return
 705:         end                
 706: !****|******************************************************************|
 707:         double precision function gfunc(x,a,b1,b2,x1,x2)
 708:         double precision x,a,b1,b2,x1,x2
 709:         
 710:         gfunc=a*exp(-b1*(x-x1)**2)+exp(-b2*(x-x2))
 711:                 
 712:         return
 713:         end                         
 714: !****|******************************************************************|
 715:         double precision function gfunc_d(x,a,b1,b2,x1,x2)
 716:         double precision x,a,b1,b2,x1,x2
 717:         
 718:         gfunc_d=-2*b1*(x-x1)*a*exp(-b1*(x-x1)**2)-b2*exp(-b2*(x-x2))
 719:                 
 720:         return
 721:         end        
 722: !****|******************************************************************|
 723:         double precision function gfunc_dd(x,a,b1,b2,x1,x2)
 724:         double precision x,a,b1,b2,x1,x2
 725:         
 726:         gfunc_dd=-2*a*b1*exp(-b1*(x-x1)**2)+4*a*b1**2*(x-x1)**2*exp(-b1*(x-x1)**2)   &
 727:                         +b2**2*exp(-b2*(x-x2))                 
 728:         return
 729:         end        
 730:                          
 731: !****|******************************************************************|
 732: !****|******************************************************************|
 733: !****|******************************************************************|
 734: 
 735: 
 736: !****|******************************************************************|
 737: !****|******************************************************************|
 738:        double precision function fpsi(x)
 739:        USE COMMONS, ONLY : ACK1, ACK2
 740:        implicit double precision (a-h,o-z) 
 741:        include 'ackland_mendelev_fe.h'
 742:        
 743:        temp = 0.d0
 744:        do i=1,npsi
 745:          temp = temp + ap(i)*Hfunc(rp(i)-x)*(rp(i)-x)**3
 746:        end do
 747:        fpsi = temp
 748:        return
 749:        end  
 750: !****|******************************************************************|
 751:        double precision function fpsi_d(x)
 752:        USE COMMONS, ONLY : ACK1, ACK2
 753:        implicit double precision (a-h,o-z) 
 754:        include 'ackland_mendelev_fe.h'
 755:        
 756:        temp = 0.d0
 757:        do i=1,npsi
 758:          temp = temp - 3.d0*ap(i)*Hfunc(rp(i)-x)*(rp(i)-x)**2
 759:        end do
 760:        fpsi_d = temp
 761:        return
 762:        end  
 763: !****|******************************************************************|
 764:        double precision function fpsi_dd(x)
 765:        USE COMMONS, ONLY : ACK1, ACK2
 766:        implicit double precision (a-h,o-z) 
 767:        include 'ackland_mendelev_fe.h'
 768:        
 769:        temp = 0.d0
 770:        do i=1,npsi
 771:          temp = temp + 6.d0*ap(i)*Hfunc(rp(i)-x)*(rp(i)-x)
 772:        end do
 773:        fpsi_dd = temp
 774:        return
 775:        end  
 776: !****|******************************************************************|
 777: !****|******************************************************************|
 778:        double precision function fvarphi (x)
 779:        USE COMMONS, ONLY : ACK1, ACK2
 780:        implicit double precision (a-h,o-z)
 781:        include 'ackland_mendelev_fe.h'
 782:         
 783:          threei   =1.d0/3.d0      
 784:         ZnFE2   = ZnFE*ZnFE
 785:        AU_TO_EV = hart
 786:         AU_TO_A = abohr
 787:          temp   = 0.d0
 788: 
 789:        rs = 0.88534d0*abohr*ZnFE**(-threei)/dsqrt(2.d0)
 790:        rx = x/rs
 791:        
 792:         if (x.lt.r1) then 
 793:           fvarphi = ZnFE2 * fphi(rx) * AU_TO_A* AU_TO_EV/x 
 794:         else if ((x.ge.r1).and.(x.lt.r2)) then 
 795:           fvarphi = dexp (  bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x )
 796:         else if (x.ge.r2) then
 797:           do i=1,nvarphi
 798:              temp = temp + af(i)*Hfunc(rf(i)-x)*(rf(i)-x)**3
 799:           end do
 800:           fvarphi = temp
 801:         end if              
 802:         
 803:        return
 804:        end 
 805: !****|******************************************************************|
 806:        double precision function fvarphi_d (x)
 807:        USE COMMONS, ONLY : ACK1, ACK2
 808:        implicit double precision (a-h,o-z)
 809:        include 'ackland_mendelev_fe.h'
 810:          threei   = 1.d0/3.d0
 811:         ZnFE2   = ZnFE*ZnFE
 812:        AU_TO_EV = hart
 813:         AU_TO_A = abohr
 814:          temp   = 0.d0
 815: 
 816:        rs = 0.88534d0*abohr*ZnFE**(-threei)/dsqrt(2.d0)
 817:        rx = x/rs
 818: 
 819:         if (x.lt.r1) then 
 820:           fvarphi_d =  - ZnFE2*fphi(rx)*AU_TO_EV*AU_TO_A/(x*x) &
 821:                  + ZnFE2 * fphi_d(rx) * AU_TO_EV*AU_TO_A / (x*rs)     
 822:         else if ((x.ge.r1).and.(x.lt.r2)) then 
 823:           fvarphi_d = (bFE1 + 2.0d0*bFE2*x + 3.0d0*bFE3*x*x) &
 824:                      *dexp (bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x)
 825:         else if (x.ge.r2) then
 826:           do i=1,nvarphi
 827:              temp = temp - 3.d0*af(i)*Hfunc(rf(i)-x)*(rf(i)-x)**2
 828:           end do
 829:           fvarphi_d = temp
 830:         end if              
 831:         
 832:        return
 833:        end
 834: !****|******************************************************************|
 835:        double precision function fvarphi_dd (x)
 836:        USE COMMONS, ONLY : ACK1, ACK2
 837:        implicit double precision (a-h,o-z)
 838:        include 'ackland_mendelev_fe.h'
 839:          threei   = 1.d0/3.d0
 840:         ZnFE2   = ZnFE*ZnFE
 841:        AU_TO_EV = hart
 842:         AU_TO_A = abohr
 843:          temp   = 0.d0
 844: 
 845:        rs = 0.88534d0*abohr*ZnFE**(-threei)/dsqrt(2.d0)
 846:        rx = x/rs
 847: 
 848:         if (x.lt.r1) then 
 849:         fvarphi_dd = 0.5d0*ZnFE2*fphi(rx)*AU_TO_EV*AU_TO_A/(x*x*x) &
 850:                      - ZnFE2*fphi_d(rx)*AU_TO_EV*AU_TO_A/(x*x*rs)  &
 851:                      - ZnFE2 * fphi_d(rx) * AU_TO_EV*AU_TO_A / (x*x*rs) &
 852:                      + ZnFE2 * fphi_dd(rx) * AU_TO_EV*AU_TO_A / (x*rs*rs)                          
 853:         else if ((x.ge.r1).and.(x.lt.r2)) then 
 854:           fvarphi_dd = (2.0d0*bFE2 + 6.0d0*bFE3*x) &
 855:                      *dexp (bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x)&
 856:                      + &
 857:                     (bFE1 + 2.0d0*bFE2*x + 3.0d0*bFE3*x*x)**2 &
 858:                      *dexp (bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x)                     
 859:         else if (x.ge.r2) then
 860:           do i=1,nvarphi
 861:              temp = temp + 6.d0*af(i)*Hfunc(rf(i)-x)*(rf(i)-x)
 862:           end do
 863:           fvarphi_dd = temp
 864:         end if              
 865:         
 866:        return
 867:        end
 868: !****|******************************************************************|
 869: !****|******************************************************************|
 870: !****|******************************************************************|
 871:        double precision function fphi(x)
 872:        double precision x
 873:   
 874:         fphi = 0.1818d0*dexp(-3.2d0*x) & 
 875:             +  0.5099d0*dexp(-0.9423d0*x) &
 876:             +  0.2802d0*dexp(-0.4029d0*x) &
 877:             +  0.02817*dexp(-0.2016*x)
 878: 
 879:        return
 880:        end
 881: !****|******************************************************************|
 882:        double precision function fphi_d(x)
 883:        double precision x
 884:   
 885:        fphi_d = -  0.58176d0*dexp(-3.2d0*x) &
 886:                 -  0.480479d0*dexp(-0.9423d0*x) & 
 887:                 -  0.112893d0*dexp(-0.4029d0*x) &
 888:                 -  0.00567907d0*dexp(-0.2016*x)
 889: 
 890:        return
 891:        end
 892: 
 893: !****|******************************************************************|
 894:        double precision function fphi_dd(x)
 895:        double precision x
 896:   
 897:        fphi_dd =    1.86163d0*dexp(-3.2d0*x) &
 898:                 +  0.452755d0*dexp(-0.9423d0*x) & 
 899:                 +  0.0454846d0*dexp(-0.4029d0*x) &
 900:                 +  0.0011449d0*dexp(-0.2016*x)
 901: 
 902:        return
 903:        end
 904: 
 905: !****|******************************************************************|
 906: !****|******************************************************************|
 907: !****|******************************************************************|
 908: !****|******************************************************************|
 909: 
 910:       SUBROUTINE BUILD_RHO_SITE(ipot,rho_site,Rn,ndir,nat_up,ndir_max)
 911:         implicit double precision (a-h,o-z)
 912:         double precision Rn(nat_up,ndir_max,3), rho_site(nat_up),    &
 913:                           rho_temp,normR,R,Rtemp(3)
 914:         integer ndir(nat_up),nat_up,ndir_max,ipot    
 915: 
 916:   
 917:        do i=1,nat_up  
 918:           rho_temp=0.0D0
 919:            do ni=1,ndir(i)
 920:             Rtemp(:)=Rn(i,ni,:)           
 921:             R=sqrt(DOT_PRODUCT(Rtemp,Rtemp)) 
 922:             rho_temp=rho_temp+rho_pot(ipot,R)            
 923:            end do
 924: 
 925:             rho_site(i)=rho_temp        
 926:         end do
 927:  
 928:        RETURN
 929:        END
 930: !****|******************************************************************|       
 931:        SUBROUTINE BUILD_V_SITE(ipot,V_site,Rn,ndir,nat_up,ndir_max)
 932:         implicit double precision (a-h,o-z)
 933:         double precision Rn(nat_up,ndir_max,3),V_site(nat_up),   &
 934:                          V_temp,normR,R,Rtemp(3)
 935:         integer ndir(nat_up),nat_up,nat_max,ndir_max,ipot
 936: 
 937:   
 938:        do i=1,nat_up  
 939:           V_temp=0.0D0
 940:            do ni=1,ndir(i)           
 941:             Rtemp(:)=Rn(i,ni,:)           
 942:             R=sqrt(DOT_PRODUCT(Rtemp,Rtemp))                      
 943:             V_temp=V_temp+Vpot(ipot,R)            
 944:            end do
 945:               V_site(i)=V_temp        
 946:         end do
 947:  
 948:        RETURN
 949:        END      
354: 950: 
355:   ELSEIF(ipot.EQ.2) THEN951:       
356:      Vpot_dd=p*(p+1)/R**2*(Ro/R)**p*fcut(ipot,R,Rc,delta)          & 
357:           -2*p/R*(Ro/R)**p*fcut_d(ipot,R,Rc,delta)        & 
358:           +(Ro/R)**p*fcut_dd(ipot,R,Rc,delta) 
359:      Vpot_dd=a_r*Vpot_dd 
360:  
361:   ELSEIF(ipot.EQ.3) THEN 
362:      Vpot_dd=p*(p+1)/R**2*(Ro/R)**p*fcut(ipot,R,Rc,delta)           & 
363:           -2*p/R*(Ro/R)**p*fcut_d(ipot,R,Rc,delta)        & 
364:           +(Ro/R)**p*fcut_dd(ipot,R,Rc,delta) 
365:      Vpot_dd=a_r*Vpot_dd 
366:  
367:   ELSEIF(ipot.EQ.4) THEN 
368:  
369:      !****|******************************************************************| 
370:      Vpot_dd=(E1*Mfunc_dd(R,R01,alpha1)+E2*Mfunc_dd(R,R02,alpha2))*fcut(ipot,R,Rc,h)       & 
371:           +    2*(E1*Mfunc_d(R,R01,alpha1)+E2*Mfunc_d(R,R02,alpha2))*fcut_d(ipot,R,Rc,h)       & 
372:           +     (E1*Mfunc(R,R01,alpha1)+E2*Mfunc(R,R02,alpha2)+dd)*fcut_dd(ipot,R,Rc,h)        & 
373:           -      12*S1*Hfunc(Rs1-R)*(Rs1-R)**2                                                        & 
374:           -      12*S2*Hfunc(Rs2-R)*(Rs2-R)**2                                                   & 
375:           -      12*S3*Hfunc(Rs3-R)*(Rs3-R)**2 
376:  
377:      !****|******************************************************************| 
378:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6)) THEN 
379:  
380:      Vpot_dd = 0.5d0*fvarphi_dd(R) 
381:  
382:   ELSE 
383:      WRITE(*,*) 'erreur de ipot' 
384:      Vpot_dd=0.0D0 
385:   ENDIF 
386:  
387:   RETURN 
388: END FUNCTION Vpot_dd 
389:  
390:  
391: !****|******************************************************************| 
392: DOUBLE PRECISION  FUNCTION fcut(ipot,R,Rc,width) 
393:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
394:   DOUBLE PRECISION R,Rc,width,test,x 
395:  
396:   fcut=0.0D0 
397:  
398:   IF(ipot.EQ.1.OR.ipot.EQ.2.OR.ipot.EQ.3) THEN 
399:      test=(R-Rc)/width 
400:      IF(test.LT.-100.d0)  THEN 
401:         fcut=1.0d0 
402:      ELSEIF(test.GT.100.d0) THEN 
403:         fcut=0.0d0 
404:      ELSE 
405:         x=test 
406:         fcut=1.0d0/(1.0d0+EXP(x)) 
407:      ENDIF 
408:   ELSEIF(ipot.EQ.4) THEN 
409:      test=(R-Rc) 
410:      IF(test.GT.0.d0)  THEN 
411:         fcut=0.0d0 
412:      ELSE 
413:         x=test/width 
414:         fcut= x**4/(1+x**4) 
415:      ENDIF 
416:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6).OR.(ipot.EQ.12).OR.(ipot.EQ.13)) THEN 
417:      fcut=1.d0 
418:   ENDIF 
419:  
420:   RETURN 
421: END FUNCTION fcut 
422:  
423: !****|******************************************************************| 
424: DOUBLE PRECISION  FUNCTION fcut_d(ipot,R,Rc,width) 
425:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
426:   DOUBLE PRECISION R,Rc,width,test,x 
427:  
428:   fcut_d=0.0D0 
429:  
430:   IF(ipot.EQ.1.OR.ipot.EQ.2.OR.ipot.EQ.3) THEN 
431:      test=(R-Rc)/width 
432:      IF(test.LT.-100.d0)  THEN 
433:         fcut_d=0.0d0 
434:      ELSEIF(test.GT.100.d0) THEN 
435:         fcut_d=0.0d0 
436:      ELSE 
437:         x=test 
438:         fcut_d=-EXP(x)/(1+EXP(x))**2 
439:         fcut_d=fcut_d/width 
440:      ENDIF 
441:   ELSEIF(ipot.EQ.4) THEN 
442:      test=(R-Rc) 
443:      IF(test.GT.0.d0)  THEN 
444:         fcut_d=0.0d0 
445:      ELSEIF(test.LE.0.d0) THEN 
446:         x=test/width 
447:         fcut_d= 4*x**3/(1+x**4)**2 
448:         fcut_d=fcut_d/width 
449:      ENDIF 
450:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6).OR.(ipot.EQ.12).OR.(ipot.EQ.13)) THEN 
451:      fcut_d=1.d0 
452:   ENDIF 
453:  
454:   RETURN 
455: END FUNCTION fcut_d 
456:  
457: !****|******************************************************************| 
458: DOUBLE PRECISION FUNCTION fcut_dd(ipot,R,Rc,width) 
459:   IMPLICIT  DOUBLE PRECISION(a-h,o-z) 
460:   DOUBLE PRECISION R,Rc,width,test,x 
461:  
462:   fcut_dd=0.0D0 
463:  
464:   IF(ipot.EQ.1.OR.ipot.EQ.2.OR.ipot.EQ.3) THEN 
465:      test=(R-Rc)/width 
466:      IF(test.LT.-100.d0)  THEN 
467:         fcut_dd=0.0d0 
468:      ELSEIF(test.GT.100.d0) THEN 
469:         fcut_dd=0.0d0 
470:      ELSE 
471:         x=test 
472:         fcut_dd= -EXP(x)/( (1.0d0+EXP(x) ) )**2 +2.0d0*EXP(2*x)/(1.0d0+EXP(x) )**3 
473:         fcut_dd=fcut_dd/width**2 
474:      ENDIF 
475:   ELSEIF(ipot.EQ.4) THEN 
476:      test=(R-Rc) 
477:      IF(test.GT.0.d0)  THEN 
478:         fcut_dd=0.0d0 
479:      ELSEIF(test.LE.0.d0) THEN 
480:         x=test/width 
481:         fcut_dd=12*x**2/(1+x**4)**2-32*x**6/(1+x**4)**3 
482:         fcut_dd=fcut_dd/width**2 
483:      ENDIF 
484:   ELSEIF ((ipot.EQ.5).OR.(ipot.EQ.6).OR.(ipot.EQ.12).OR.(ipot.EQ.13)) THEN 
485:      fcut_dd=1.d0 
486:  
487:   ENDIF 
488:  
489:   RETURN 
490: END FUNCTION fcut_dd 
491:  
492: !****|******************************************************************| 
493: DOUBLE PRECISION FUNCTION Fembed(ipot,x) 
494:   USE COMMONS, ONLY : ACK1, ACK2 
495:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
496:   DOUBLE PRECISION x 
497:   INTEGER ipot 
498:   DOUBLE PRECISION Rc 
499:   COMMON /param_cut_off/Rc 
500:   INCLUDE 'ackland_sma.h' 
501:   INCLUDE 'ackland_mishin_cu.h' 
502:   INCLUDE 'ackland_mendelev_fe.h' 
503:  
504:   Fembed=0.0D0 
505:  
506:   IF(ipot.EQ.1.OR.ipot.EQ.2.OR.ipot.EQ.3) THEN 
507:      Fembed=x**(alpha) 
508:      Fembed=-b_a*Fembed 
509:  
510:   ELSEIF(ipot.EQ.4) THEN 
511:  
512:      IF(x.LE.1) THEN 
513:         Fembed=F0+0.5*F2*(x-1)**2+q1*(x-1)**3+q2*(x-1)**4 +q3*(x-1)**5 +q4*(x-1)**6 
514:  
515:      ELSEIF(x.GT.1) THEN 
516:         Fembed=(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)/(1+qq2*(x-1)**3) 
517:      ENDIF 
518:  
519:   ELSEIF(ipot.EQ.5) THEN 
520:  
521:      Fembed = -dsqrt(x) + aphi*x*x + aphi2*x*x*x*x 
522:  
523:   ELSEIF(ipot.EQ.6) THEN 
524:      IF (x.LT.1.d0) THEN 
525:         Fembed = -aphi*dsqrt(x) - aphi2*(1-dsqrt(x))*dlog(2.d0-x)/dlog(2.d0) 
526:      ELSE 
527:         Fembed = -aphi*dsqrt(x) 
528:      ENDIF 
529:  
530:   ELSEIF(ipot.EQ.12) THEN 
531:      Fembed = -     dsqrt(x)                                         & 
532:           - 1.9162462126235d0*1.d-7*(x-60.d0)**4*Hfunc(x-60.d0) & 
533:           + 4.6418727035037d0*1.d-7*(x-70.d0)**4*Hfunc(x-70.d0) & 
534:           + 6.6448294272955d0*1.d-7*(x-80.d0)**4*Hfunc(x-80.d0) & 
535:           - 2.0680252960229d0*1.d-6*(x-85.d0)**4*Hfunc(x-85.d0) & 
536:           + 1.1387131464983d0*1.d-6*(x-90.d0)**4*Hfunc(x-90.d0) 
537:   ELSEIF(ipot.EQ.13) THEN 
538:      Fembed = - dsqrt(x)                                         & 
539:           + 3.2283012597866d0*1.d-7*(x-60.d0)**4*Hfunc(x-60.d0) & 
540:           - 1.1552813894483d0*1.d-6*(x-70.d0)**4*Hfunc(x-70.d0) & 
541:           + 2.3747280268355d0*1.d-6*(x-80.d0)**4*Hfunc(x-80.d0) & 
542:           - 2.0379550826523d0*1.d-6*(x-85.d0)**4*Hfunc(x-85.d0) & 
543:           + 4.9758343293936d0*1.d-7*(x-90.d0)**4*Hfunc(x-90.d0) 
544:  
545:  
546:   END IF 
547:   RETURN 
548: END FUNCTION Fembed 
549: !****|******************************************************************| 
550: DOUBLE PRECISION FUNCTION Fembed_d(ipot,x) 
551:   USE COMMONS, ONLY : ACK1, ACK2 
552:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
553:   DOUBLE PRECISION x 
554:   DOUBLE PRECISION Rc 
555:   COMMON /param_cut_off/Rc 
556:   INCLUDE 'ackland_sma.h' 
557:   INCLUDE 'ackland_mishin_cu.h' 
558:   INCLUDE 'ackland_mendelev_fe.h' 
559:  
560:   Fembed_d=0.0D0 
561:  
562:   IF(ipot.EQ.1.OR.ipot.EQ.2.OR.ipot.EQ.3) THEN 
563:      Fembed_d=alpha*x**(alpha-1.0d0) 
564:      Fembed_d=-b_a*Fembed_d 
565:   ELSEIF(ipot.EQ.4) THEN 
566:      IF(x.LE.1) THEN 
567:         Fembed_d=F2*(x-1)+3*q1*(x-1)**2+4*q2*(x-1)**3+5*q3*(x-1)**4+6*q4*(x-1)**5 
568:  
569:      ELSEIF(x.GT.1) THEN 
570:         Fembed_d=(F2*(x-1)+3*q1*(x-1)**2+4*qq1*(x-1)**3)/(1+qq2*(x-1)**3)       & 
571:              -3*qq2*(x-1)**2*(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)  & 
572:              /(1+qq2*(x-1)**3)**2 
573:  
574:  
575:      ENDIF 
576:  
577:   ELSEIF (ipot.EQ.5) THEN 
578:  
579:      Fembed_d = -0.5d0/dsqrt(x) + 2.d0*aphi*x  & 
580:           + 4.d0*aphi2*x*x*x 
581:  
582:  
583:   ELSEIF(ipot.EQ.6) THEN 
584:      IF (x.LT.1.d0) THEN 
585:         Fembed_d = -aphi*0.5d0/dsqrt(x)                   & 
586:              + aphi2*(0.5d0*dlog(2.d0-x)/dsqrt(x)  & 
587:              + (1-dsqrt(x))/(2.d0-x))/dlog(2.d0) 
588:      ELSE 
589:         Fembed_d = -0.5d0*aphi/dsqrt(x) 
590:      END IF 
591:  
592:   ELSEIF(ipot.EQ.12) THEN 
593:      Fembed_d = - 0.5d0 / dsqrt(x)                                 & 
594:           + 4.d0*(                                           & 
595:           - 1.9162462126235d0*1.d-7*(x-60.d0)**3*Hfunc(x-60.d0) & 
596:           + 4.6418727035037d0*1.d-7*(x-70.d0)**3*Hfunc(x-70.d0) & 
597:           + 6.6448294272955d0*1.d-7*(x-80.d0)**3*Hfunc(x-80.d0) & 
598:           - 2.0680252960229d0*1.d-6*(x-85.d0)**3*Hfunc(x-85.d0) & 
599:           + 1.1387131464983d0*1.d-6*(x-90.d0)**3*Hfunc(x-90.d0) & 
600:           ) 
601:   ELSEIF(ipot.EQ.13) THEN 
602:      Fembed_d = - 0.5d0 / dsqrt(x)                                 & 
603:           + 4.d0*(                                           & 
604:           + 3.2283012597866d0*1.d-7*(x-60.d0)**3*Hfunc(x-60.d0) & 
605:           - 1.1552813894483d0*1.d-6*(x-70.d0)**3*Hfunc(x-70.d0) & 
606:           + 2.3747280268355d0*1.d-6*(x-80.d0)**3*Hfunc(x-80.d0) & 
607:           - 2.0379550826523d0*1.d-6*(x-85.d0)**3*Hfunc(x-85.d0) & 
608:           + 4.9758343293936d0*1.d-7*(x-90.d0)**3*Hfunc(x-90.d0) & 
609:           ) 
610:  
611:  
612:   ENDIF 
613:  
614:   RETURN 
615: END FUNCTION Fembed_d 
616:  
617: !****|******************************************************************| 
618: DOUBLE PRECISION FUNCTION Fembed_dd(ipot,x) 
619:   USE COMMONS, ONLY : ACK1, ACK2 
620:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
621:   INCLUDE 'ackland_sma.h' 
622:   INCLUDE 'ackland_mishin_cu.h' 
623:   DOUBLE PRECISION Rc 
624:   COMMON /param_cut_off/Rc 
625:   DOUBLE PRECISION x,denom 
626:   INCLUDE 'ackland_mendelev_fe.h' 
627:  
628:   Fembed_dd=0.0D0 
629:  
630:   IF(ipot.EQ.1.OR.ipot.EQ.2.OR.ipot.EQ.3) THEN 
631:      Fembed_dd=alpha*(alpha-1)*x**(alpha-2.0d0) 
632:      Fembed_dd=-b_a*Fembed_dd 
633:   ELSEIF(ipot.EQ.4) THEN 
634:      IF(x<1) THEN 
635:         Fembed_dd=F2+6*q1*(x-1)+12*q2*(x-1)**2+20*q3*(x-1)**3+30*q4*(x-1)**4 
636:  
637:      ELSEIF(x>1) THEN 
638:         denom=1+qq2*(x-1)**3 
639:  
640:         Fembed_dd=(F2+6*q1*(x-1)+12*qq1*(x-1)**2)/denom                              & 
641:              -6*qq2*(x-1)**2*(F2*(x-1)+3*q1*(x-1)**2+4*qq1*(x-1)**3)/denom**2             & 
642:              -6*qq2*(x-1)*(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)/denom**2          & 
643:              +18*qq2**2*(x-1)**4*(F0+0.5*F2*(x-1)**2+q1*(x-1)**3+qq1*(x-1)**4)/denom**3 
644:  
645:      ENDIF 
646:  
647:   ELSEIF (ipot.EQ.5) THEN 
648:      Fembed_dd = 0.25d0/dsqrt(x**3) + 2.d0*aphi  & 
649:           + 12.d0*aphi2*x*x 
650:  
651:   ELSEIF (ipot.EQ.6) THEN 
652:      IF (x.LT.1.d0) THEN 
653:         Fembed_dd = 0.25d0*aphi/dsqrt(x**3)  + aphi2/dlog(2.d0)*   & 
654:              (                                              & 
655:              -0.25d0*dlog(2.d0-x)/dsqrt(x**3)             & 
656:              -1.0/((2.d0-x)*dsqrt(x))                     & 
657:              +(1.0-dsqrt(x))/(2.d0-x)**2                  & 
658:              ) 
659:      ELSE 
660:         Fembed_dd = 0.25d0*aphi/dsqrt(x**3) 
661:      END IF 
662:  
663:   ENDIF 
664:  
665:   RETURN 
666: END FUNCTION Fembed_dd 
667: !****|******************************************************************| 
668: DOUBLE PRECISION FUNCTION delta_dirac(i,j) 
669:   INTEGER i,j 
670:  
671:   delta_dirac=0.0d0 
672:   IF(i.EQ.j) THEN 
673:      delta_dirac=1.0d0 
674:   ENDIF 
675:   RETURN 
676: END FUNCTION delta_dirac 
677: !****|******************************************************************| 
678: DOUBLE PRECISION FUNCTION Hfunc(x) 
679:   DOUBLE PRECISION x 
680:  
681:   IF(x.LT.0.0) THEN 
682:      Hfunc=0.0d0 
683:   ELSE 
684:      Hfunc=1.0d0 
685:   ENDIF 
686:   RETURN 
687: END FUNCTION Hfunc 
688: !****|******************************************************************| 
689: DOUBLE PRECISION FUNCTION Hfunc0(x) 
690:   DOUBLE PRECISION x 
691:  
692:   IF(x.LE.0.0) THEN 
693:      Hfunc0=0.0d0 
694:   ELSE 
695:      Hfunc0=1.0d0 
696:   ENDIF 
697:   RETURN 
698: END FUNCTION Hfunc0 
699: !****|******************************************************************| 
700: DOUBLE PRECISION FUNCTION Mfunc(R,R0,alpha) 
701:   DOUBLE PRECISION R,R0,alpha 
702:  
703:   Mfunc=EXP(-2*alpha*(R-R0))-2*EXP(-alpha*(R-R0)) 
704:  
705:   RETURN 
706: END FUNCTION Mfunc 
707: !****|******************************************************************| 
708: DOUBLE PRECISION FUNCTION Mfunc_d(R,R0,alpha) 
709:   DOUBLE PRECISION R,R0,alpha 
710:  
711:   Mfunc_d=-2*alpha*EXP(-2*alpha*(R-R0))+2*alpha*EXP(-alpha*(R-R0)) 
712:  
713:   RETURN 
714: END FUNCTION Mfunc_d 
715: !****|******************************************************************| 
716: DOUBLE PRECISION FUNCTION Mfunc_dd(R,R0,alpha) 
717:   DOUBLE PRECISION R,R0,alpha 
718:  
719:   Mfunc_dd=+4*alpha**2*EXP(-2*alpha*(R-R0))-2*alpha**2*EXP(-alpha*(R-R0)) 
720:  
721:   RETURN 
722: END FUNCTION Mfunc_dd 
723: !****|******************************************************************| 
724: DOUBLE PRECISION FUNCTION gfunc(x,a,b1,b2,x1,x2) 
725:   DOUBLE PRECISION x,a,b1,b2,x1,x2 
726:  
727:   gfunc=a*EXP(-b1*(x-x1)**2)+EXP(-b2*(x-x2)) 
728:  
729:   RETURN 
730: END FUNCTION gfunc 
731: !****|******************************************************************| 
732: DOUBLE PRECISION FUNCTION gfunc_d(x,a,b1,b2,x1,x2) 
733:   DOUBLE PRECISION x,a,b1,b2,x1,x2 
734:  
735:   gfunc_d=-2*b1*(x-x1)*a*EXP(-b1*(x-x1)**2)-b2*EXP(-b2*(x-x2)) 
736:  
737:   RETURN 
738: END FUNCTION gfunc_d 
739: !****|******************************************************************| 
740: DOUBLE PRECISION FUNCTION gfunc_dd(x,a,b1,b2,x1,x2) 
741:   DOUBLE PRECISION x,a,b1,b2,x1,x2 
742:  
743:   gfunc_dd=-2*a*b1*EXP(-b1*(x-x1)**2)+4*a*b1**2*(x-x1)**2*EXP(-b1*(x-x1)**2)   & 
744:        +b2**2*EXP(-b2*(x-x2)) 
745:   RETURN 
746: END FUNCTION gfunc_dd 
747:  
748: !****|******************************************************************| 
749: !****|******************************************************************| 
750: !****|******************************************************************| 
751:  
752:  
753: !****|******************************************************************| 
754: !****|******************************************************************| 
755: DOUBLE PRECISION FUNCTION fpsi(x) 
756:   USE COMMONS, ONLY : ACK1, ACK2 
757:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
758:   INCLUDE 'ackland_mendelev_fe.h' 
759:  
760:   temp = 0.d0 
761:   DO i=1,npsi 
762:      temp = temp + ap(i)*Hfunc(rp(i)-x)*(rp(i)-x)**3 
763:   END DO 
764:   fpsi = temp 
765:   RETURN 
766: END FUNCTION fpsi 
767: !****|******************************************************************| 
768: DOUBLE PRECISION FUNCTION fpsi_d(x) 
769:   USE COMMONS, ONLY : ACK1, ACK2 
770:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
771:   INCLUDE 'ackland_mendelev_fe.h' 
772:  
773:   temp = 0.d0 
774:   DO i=1,npsi 
775:      temp = temp - 3.d0*ap(i)*Hfunc(rp(i)-x)*(rp(i)-x)**2 
776:   END DO 
777:   fpsi_d = temp 
778:   RETURN 
779: END FUNCTION fpsi_d 
780: !****|******************************************************************| 
781: DOUBLE PRECISION FUNCTION fpsi_dd(x) 
782:   USE COMMONS, ONLY : ACK1, ACK2 
783:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
784:   INCLUDE 'ackland_mendelev_fe.h' 
785:  
786:   temp = 0.d0 
787:   DO i=1,npsi 
788:      temp = temp + 6.d0*ap(i)*Hfunc(rp(i)-x)*(rp(i)-x) 
789:   END DO 
790:   fpsi_dd = temp 
791:   RETURN 
792: END FUNCTION fpsi_dd 
793: !****|******************************************************************| 
794: !****|******************************************************************| 
795: DOUBLE PRECISION FUNCTION fvarphi (x) 
796:   USE COMMONS, ONLY : ACK1, ACK2 
797:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
798:   INCLUDE 'ackland_mendelev_fe.h' 
799:  
800:   threei   =1.d0/3.d0 
801:   ZnFE2   = ZnFE*ZnFE 
802:   AU_TO_EV = hart 
803:   AU_TO_A = abohr 
804:   temp   = 0.d0 
805:  
806:   rs = 0.88534d0*abohr*ZnFE**(-threei)/dsqrt(2.d0) 
807:   rx = x/rs 
808:  
809:   IF (x.LT.r1) THEN 
810:      fvarphi = ZnFE2 * fphi(rx) * AU_TO_A* AU_TO_EV/x 
811:   ELSE IF ((x.GE.r1).AND.(x.LT.r2)) THEN 
812:      fvarphi = dexp (  bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x ) 
813:   ELSE 
814:      DO i=1,nvarphi 
815:         temp = temp + af(i)*Hfunc(rf(i)-x)*(rf(i)-x)**3 
816:      END DO 
817:      fvarphi = temp 
818:   END IF 
819:  
820:   RETURN 
821: END FUNCTION fvarphi 
822: !****|******************************************************************| 
823: DOUBLE PRECISION FUNCTION fvarphi_d (x) 
824:   USE COMMONS, ONLY : ACK1, ACK2 
825:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
826:   INCLUDE 'ackland_mendelev_fe.h' 
827:   threei   = 1.d0/3.d0 
828:   ZnFE2   = ZnFE*ZnFE 
829:   AU_TO_EV = hart 
830:   AU_TO_A = abohr 
831:   temp   = 0.d0 
832:  
833:   rs = 0.88534d0*abohr*ZnFE**(-threei)/dsqrt(2.d0) 
834:   rx = x/rs 
835:  
836:   IF (x.LT.r1) THEN 
837:      fvarphi_d =  - ZnFE2*fphi(rx)*AU_TO_EV*AU_TO_A/(x*x) & 
838:           + ZnFE2 * fphi_d(rx) * AU_TO_EV*AU_TO_A / (x*rs) 
839:   ELSE IF ((x.GE.r1).AND.(x.LT.r2)) THEN 
840:      fvarphi_d = (bFE1 + 2.0d0*bFE2*x + 3.0d0*bFE3*x*x) & 
841:           *dexp (bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x) 
842:   ELSE 
843:      DO i=1,nvarphi 
844:         temp = temp - 3.d0*af(i)*Hfunc(rf(i)-x)*(rf(i)-x)**2 
845:      END DO 
846:      fvarphi_d = temp 
847:   END IF 
848:  
849:   RETURN 
850: END FUNCTION fvarphi_d 
851: !****|******************************************************************| 
852: DOUBLE PRECISION FUNCTION fvarphi_dd (x) 
853:   USE COMMONS, ONLY : ACK1, ACK2 
854:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
855:   INCLUDE 'ackland_mendelev_fe.h' 
856:   threei   = 1.d0/3.d0 
857:   ZnFE2   = ZnFE*ZnFE 
858:   AU_TO_EV = hart 
859:   AU_TO_A = abohr 
860:   temp   = 0.d0 
861:  
862:   rs = 0.88534d0*abohr*ZnFE**(-threei)/dsqrt(2.d0) 
863:   rx = x/rs 
864:  
865:   IF (x.LT.r1) THEN 
866:      fvarphi_dd = 0.5d0*ZnFE2*fphi(rx)*AU_TO_EV*AU_TO_A/(x*x*x) & 
867:           - ZnFE2*fphi_d(rx)*AU_TO_EV*AU_TO_A/(x*x*rs)  & 
868:           - ZnFE2 * fphi_d(rx) * AU_TO_EV*AU_TO_A / (x*x*rs) & 
869:           + ZnFE2 * fphi_dd(rx) * AU_TO_EV*AU_TO_A / (x*rs*rs) 
870:   ELSE IF ((x.GE.r1).AND.(x.LT.r2)) THEN 
871:      fvarphi_dd = (2.0d0*bFE2 + 6.0d0*bFE3*x) & 
872:           *dexp (bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x)& 
873:           + & 
874:           (bFE1 + 2.0d0*bFE2*x + 3.0d0*bFE3*x*x)**2 & 
875:           *dexp (bFE0 + bFE1*x + bFE2*x*x + bFE3*x*x*x) 
876:   ELSE 
877:      DO i=1,nvarphi 
878:         temp = temp + 6.d0*af(i)*Hfunc(rf(i)-x)*(rf(i)-x) 
879:      END DO 
880:      fvarphi_dd = temp 
881:   END IF 
882:  
883:   RETURN 
884: END FUNCTION fvarphi_dd 
885: !****|******************************************************************| 
886: !****|******************************************************************| 
887: !****|******************************************************************| 
888: DOUBLE PRECISION FUNCTION fphi(x) 
889:   DOUBLE PRECISION x 
890:  
891:   fphi = 0.1818d0*dexp(-3.2d0*x) & 
892:        +  0.5099d0*dexp(-0.9423d0*x) & 
893:        +  0.2802d0*dexp(-0.4029d0*x) & 
894:        +  0.02817*dexp(-0.2016*x) 
895:  
896:   RETURN 
897: END FUNCTION fphi 
898: !****|******************************************************************| 
899: DOUBLE PRECISION FUNCTION fphi_d(x) 
900:   DOUBLE PRECISION x 
901:  
902:   fphi_d = -  0.58176d0*dexp(-3.2d0*x) & 
903:        -  0.480479d0*dexp(-0.9423d0*x) & 
904:        -  0.112893d0*dexp(-0.4029d0*x) & 
905:        -  0.00567907d0*dexp(-0.2016*x) 
906:  
907:   RETURN 
908: END FUNCTION fphi_d 
909:  
910: !****|******************************************************************| 
911: DOUBLE PRECISION FUNCTION fphi_dd(x) 
912:   DOUBLE PRECISION x 
913:  
914:   fphi_dd =    1.86163d0*dexp(-3.2d0*x) & 
915:        +  0.452755d0*dexp(-0.9423d0*x) & 
916:        +  0.0454846d0*dexp(-0.4029d0*x) & 
917:        +  0.0011449d0*dexp(-0.2016*x) 
918:  
919:   RETURN 
920: END FUNCTION fphi_dd 
921:  
922: !****|******************************************************************| 
923: !****|******************************************************************| 
924: !****|******************************************************************| 
925: !****|******************************************************************| 
926:  
927: SUBROUTINE BUILD_RHO_SITE(ipot,rho_site,Rn,ndir,nat_up,ndir_max) 
928:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
929:   DOUBLE PRECISION Rn(nat_up,ndir_max,3), rho_site(nat_up),    & 
930:        rho_temp,normR,R,Rtemp(3) 
931:   INTEGER ndir(nat_up),nat_up,ndir_max,ipot 
932:  
933:  
934:   DO i=1,nat_up 
935:      rho_temp=0.0D0 
936:      DO ni=1,ndir(i) 
937:         Rtemp(:)=Rn(i,ni,:) 
938:         R=SQRT(DOT_PRODUCT(Rtemp,Rtemp)) 
939:         rho_temp=rho_temp+rho_pot(ipot,R) 
940:      END DO 
941:  
942:      rho_site(i)=rho_temp 
943:   END DO 
944:  
945:   RETURN 
946: END SUBROUTINE BUILD_RHO_SITE 
947: !****|******************************************************************| 
948: SUBROUTINE BUILD_V_SITE(ipot,V_site,Rn,ndir,nat_up,ndir_max) 
949:   IMPLICIT DOUBLE PRECISION (a-h,o-z) 
950:   DOUBLE PRECISION Rn(nat_up,ndir_max,3),V_site(nat_up),   & 
951:        V_temp,normR,R,Rtemp(3) 
952:   INTEGER ndir(nat_up),nat_up,nat_max,ndir_max,ipot 
953:  
954:  
955:   DO i=1,nat_up 
956:      V_temp=0.0D0 
957:      DO ni=1,ndir(i) 
958:         Rtemp(:)=Rn(i,ni,:) 
959:         R=SQRT(DOT_PRODUCT(Rtemp,Rtemp)) 
960:         V_temp=V_temp+Vpot(ipot,R) 
961:      END DO 
962:      V_site(i)=V_temp 
963:   END DO 
964: 952: 
965:   RETURN 
966: END SUBROUTINE BUILD_V_SITE 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0