UseTheForce/DarkSide/switcheroo.F90

!!
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!
!! SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
!! research, please cite the appropriate papers when you publish your
!! work.  Good starting points are:
!!                                                                      
!! [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
!! [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
!! [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
!! [4]  Vardeman & Gezelter, in progress (2009).
!!
!!$
module switcheroo

  use definitions
  use interpolation
  use status

  implicit none
  PRIVATE

#define __FORTRAN90
#include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
  
  !! number of points for the spline approximations
  INTEGER, PARAMETER :: np = 150

  real ( kind = dp ), save :: rin
  real ( kind = dp ), save :: rout
  real ( kind = dp ), save :: rin2
  real ( kind = dp ), save :: rout2

  logical, save :: haveSplines = .false.
  logical, save :: switchIsCubic = .true.
  integer, save :: functionType = CUBIC

  ! spline structures
  type(cubicSpline), save :: r2Spline
  type(cubicSpline), save :: switchSpline

  public::set_switch_type
  public::set_switch
  public::get_switch
  public::delete_switch

contains

  subroutine set_switch(rinner, router)

    real ( kind = dp ), intent(in):: rinner, router
    real ( kind = dp ), dimension(np) :: xvals, yvals
    real ( kind = dp ), dimension(2) :: rCubVals, sCubVals
    real ( kind = dp ) :: rval, rval2, rval3, rval4, rval5
    real ( kind = dp ) :: rvaldi, rvaldi2, rvaldi3, rvaldi4, rvaldi5
    real ( kind = dp ) :: c0, c3, c4, c5, dx, r, r2
    integer :: i

    if (router .lt. rinner) then
       call handleError("set_switch", "router is less than rinner")
       return
    endif

    if ((router .lt. 0.0_dp) .or. (rinner .lt. 0.0_dp))  then
       call handleError("set_switch", "one of the switches is negative!")
       return
    endif

    rin = rinner
    rout = router
    rin2 = rinner * rinner
    rout2 = router * router

    if ((router-rinner) .lt. 1e-8)  then
       ! no reason to set up splines if the switching region is tiny
       return
    endif

    dx = (rout2-rin2) / dble(np-1)
    
    do i = 1, np
       r2 = rin2 + dble(i-1)*dx
       xvals(i) = r2
       yvals(i) = sqrt(r2)
    enddo

    call newSpline(r2spline, xvals, yvals, .true.)

    if (functionType .eq. FIFTH_ORDER_POLY) then
       c0 = 1.0_dp
       c3 = -10.0_dp
       c4 = 15.0_dp
       c5 = -6.0_dp

       dx = (rout-rin) / dble(np-1)
    
       do i = 1, np
          r = rin + dble(i-1)*dx
          xvals(i) = r

          rval = ( r - rin )
          rval2 = rval*rval
          rval3 = rval2*rval
          rval4 = rval2*rval2
          rval5 = rval3*rval2
          rvaldi = 1.0_dp/( rout - rin )
          rvaldi2 = rvaldi*rvaldi
          rvaldi3 = rvaldi2*rvaldi
          rvaldi4 = rvaldi2*rvaldi2
          rvaldi5 = rvaldi3*rvaldi2
          yvals(i)= c0 + c3*rval3*rvaldi3 + c4*rval4*rvaldi4 + c5*rval5*rvaldi5
       enddo
       
       call newSpline(switchSpline, xvals, yvals, .true.)
       
       switchIsCubic = .false.
    else
       rCubVals(1) = rin
       rCubVals(2) = rout
       sCubVals(1) = 1.0_dp
       sCubVals(2) = 0.0_dp      
       call newSpline(switchSpline, rCubVals, sCubVals, .true.)
    endif
    
    haveSplines = .true.
    return
  end subroutine set_switch

  subroutine set_switch_type(functionForm)
    integer, intent(in) :: functionForm    
    functionType = functionForm

    if ((functionType.eq.FIFTH_ORDER_POLY).or.(functionType.eq.CUBIC)) then
       if (haveSplines) then
          call delete_switch()
          call set_switch(rin, rout)
       endif
    else
       call handleError("set_switch_type", &
            "Unknown type of switching function!")
       return      
    endif
  end subroutine set_switch_type
  
  subroutine delete_switch()
    call deleteSpline(switchSpline)
    call deleteSpline(r2spline)
    return
  end subroutine delete_switch
  
  subroutine get_switch(r2, sw, dswdr, r, in_switching_region)

    real( kind = dp ), intent(in) :: r2
    real( kind = dp ), intent(inout) :: sw, dswdr, r
    logical, intent(inout) :: in_switching_region
    integer :: j
    real ( kind = dp ) :: a, b, c, d, dx

    sw = 1.0_dp
    dswdr = 0.0_dp
    in_switching_region = .false.

    if (r2.gt.rin2) then
       if (r2.gt.rout2) then

          sw = 0.0_dp
          dswdr = 0.0_dp
          return
          
       else         
          
          call lookupUniformSpline(r2Spline, r2, r)
          if (switchIsCubic) then
             ! super zippy automated use of precomputed spline coefficients
             dx = r - rin
             sw = switchSpline%y(1) + dx*(dx*(switchSpline%c(1) + &
                  dx*switchSpline%d(1)))
             dswdr = dx*(2.0_dp * switchSpline%c(1) + &
                  3.0_dp * dx * switchSpline%d(1))
          else
             call lookupUniformSpline1d(switchSpline, r, sw, dswdr)
          endif
          
          in_switching_region = .true.
          
          return          
       endif
    else
       return
    endif
    
  end subroutine get_switch

end module switcheroo
Revision:	1465
Committed:	Fri Jul 9 23:08:25 2010 UTC (14 years, 11 months ago) by chuckv
File size:	6796 byte(s)
Log Message:	Creating busticated version of OpenMD
#	User	Rev	Content
1	gezelter	246	!!
2			!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9	gezelter	1390	!! 1. Redistributions of source code must retain the above copyright
10	gezelter	246	!! notice, this list of conditions and the following disclaimer.
11			!!
12	gezelter	1390	!! 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	!! notice, this list of conditions and the following disclaimer in the
14			!! documentation and/or other materials provided with the
15			!! distribution.
16			!!
17			!! This software is provided "AS IS," without a warranty of any
18			!! kind. All express or implied conditions, representations and
19			!! warranties, including any implied warranty of merchantability,
20			!! fitness for a particular purpose or non-infringement, are hereby
21			!! excluded. The University of Notre Dame and its licensors shall not
22			!! be liable for any damages suffered by licensee as a result of
23			!! using, modifying or distributing the software or its
24			!! derivatives. In no event will the University of Notre Dame or its
25			!! licensors be liable for any lost revenue, profit or data, or for
26			!! direct, indirect, special, consequential, incidental or punitive
27			!! damages, however caused and regardless of the theory of liability,
28			!! arising out of the use of or inability to use software, even if the
29			!! University of Notre Dame has been advised of the possibility of
30			!! such damages.
31			!!
32	gezelter	1390	!! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			!! research, please cite the appropriate papers when you publish your
34			!! work. Good starting points are:
35			!!
36			!! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			!! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			!! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39			!! [4] Vardeman & Gezelter, in progress (2009).
40			!!
41	gezelter	939	!!$
42	gezelter	115	module switcheroo
43
44			use definitions
45	chrisfen	937	use interpolation
46	gezelter	939	use status
47	gezelter	115
48			implicit none
49			PRIVATE
50
51			#define __FORTRAN90
52	chrisfen	725	#include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
53	gezelter	939
54			!! number of points for the spline approximations
55			INTEGER, PARAMETER :: np = 150
56	gezelter	115
57	gezelter	939	real ( kind = dp ), save :: rin
58			real ( kind = dp ), save :: rout
59			real ( kind = dp ), save :: rin2
60			real ( kind = dp ), save :: rout2
61	gezelter	115
62	gezelter	939	logical, save :: haveSplines = .false.
63			logical, save :: switchIsCubic = .true.
64	chrisfen	725	integer, save :: functionType = CUBIC
65	gezelter	115
66	gezelter	939	! spline structures
67			type(cubicSpline), save :: r2Spline
68			type(cubicSpline), save :: switchSpline
69	gezelter	115
70	gezelter	939	public::set_switch_type
71	gezelter	115	public::set_switch
72			public::get_switch
73	gezelter	939	public::delete_switch
74	gezelter	115
75			contains
76
77	gezelter	939	subroutine set_switch(rinner, router)
78	gezelter	938
79	gezelter	115	real ( kind = dp ), intent(in):: rinner, router
80	gezelter	939	real ( kind = dp ), dimension(np) :: xvals, yvals
81			real ( kind = dp ), dimension(2) :: rCubVals, sCubVals
82			real ( kind = dp ) :: rval, rval2, rval3, rval4, rval5
83			real ( kind = dp ) :: rvaldi, rvaldi2, rvaldi3, rvaldi4, rvaldi5
84			real ( kind = dp ) :: c0, c3, c4, c5, dx, r, r2
85	chrisfen	937	integer :: i
86	gezelter	115
87			if (router .lt. rinner) then
88	gezelter	939	call handleError("set_switch", "router is less than rinner")
89	gezelter	115	return
90			endif
91
92	gezelter	960	if ((router .lt. 0.0_dp) .or. (rinner .lt. 0.0_dp)) then
93	gezelter	939	call handleError("set_switch", "one of the switches is negative!")
94	gezelter	115	return
95			endif
96	gezelter	507
97	gezelter	939	rin = rinner
98			rout = router
99			rin2 = rinner * rinner
100			rout2 = router * router
101	gezelter	115
102	gezelter	939	if ((router-rinner) .lt. 1e-8) then
103			! no reason to set up splines if the switching region is tiny
104			return
105	chrisfen	937	endif
106	gezelter	939
107			dx = (rout2-rin2) / dble(np-1)
108	chrisfen	937
109	gezelter	939	do i = 1, np
110			r2 = rin2 + dble(i-1)*dx
111			xvals(i) = r2
112	gezelter	960	yvals(i) = sqrt(r2)
113	gezelter	939	enddo
114	gezelter	115
115	gezelter	939	call newSpline(r2spline, xvals, yvals, .true.)
116	chrisfen	725
117			if (functionType .eq. FIFTH_ORDER_POLY) then
118	gezelter	960	c0 = 1.0_dp
119			c3 = -10.0_dp
120			c4 = 15.0_dp
121			c5 = -6.0_dp
122	gezelter	939
123			dx = (rout-rin) / dble(np-1)
124
125			do i = 1, np
126			r = rin + dble(i-1)*dx
127			xvals(i) = r
128
129			rval = ( r - rin )
130			rval2 = rval*rval
131			rval3 = rval2*rval
132			rval4 = rval2*rval2
133			rval5 = rval3*rval2
134	gezelter	960	rvaldi = 1.0_dp/( rout - rin )
135	gezelter	939	rvaldi2 = rvaldi*rvaldi
136			rvaldi3 = rvaldi2*rvaldi
137			rvaldi4 = rvaldi2*rvaldi2
138			rvaldi5 = rvaldi3*rvaldi2
139			yvals(i)= c0 + c3rval3rvaldi3 + c4rval4rvaldi4 + c5rval5rvaldi5
140			enddo
141
142			call newSpline(switchSpline, xvals, yvals, .true.)
143
144			switchIsCubic = .false.
145			else
146			rCubVals(1) = rin
147			rCubVals(2) = rout
148	gezelter	960	sCubVals(1) = 1.0_dp
149			sCubVals(2) = 0.0_dp
150	gezelter	939	call newSpline(switchSpline, rCubVals, sCubVals, .true.)
151	chrisfen	725	endif
152	gezelter	939
153			haveSplines = .true.
154			return
155			end subroutine set_switch
156	chrisfen	725
157	gezelter	939	subroutine set_switch_type(functionForm)
158			integer, intent(in) :: functionForm
159			functionType = functionForm
160	gezelter	115
161	gezelter	939	if ((functionType.eq.FIFTH_ORDER_POLY).or.(functionType.eq.CUBIC)) then
162			if (haveSplines) then
163			call delete_switch()
164			call set_switch(rin, rout)
165			endif
166			else
167			call handleError("set_switch_type", &
168			"Unknown type of switching function!")
169			return
170			endif
171			end subroutine set_switch_type
172
173			subroutine delete_switch()
174			call deleteSpline(switchSpline)
175			call deleteSpline(r2spline)
176			return
177			end subroutine delete_switch
178
179			subroutine get_switch(r2, sw, dswdr, r, in_switching_region)
180
181	gezelter	115	real( kind = dp ), intent(in) :: r2
182			real( kind = dp ), intent(inout) :: sw, dswdr, r
183			logical, intent(inout) :: in_switching_region
184	gezelter	938	integer :: j
185	gezelter	939	real ( kind = dp ) :: a, b, c, d, dx
186	gezelter	115
187	gezelter	960	sw = 1.0_dp
188			dswdr = 0.0_dp
189	gezelter	115	in_switching_region = .false.
190
191	chrisfen	940	if (r2.gt.rin2) then
192			if (r2.gt.rout2) then
193	gezelter	115
194	gezelter	960	sw = 0.0_dp
195			dswdr = 0.0_dp
196	gezelter	115	return
197	chrisfen	940
198	gezelter	939	else
199
200			call lookupUniformSpline(r2Spline, r2, r)
201			if (switchIsCubic) then
202			! super zippy automated use of precomputed spline coefficients
203			dx = r - rin
204			sw = switchSpline%y(1) + dx(dx(switchSpline%c(1) + &
205			dx*switchSpline%d(1)))
206	gezelter	960	dswdr = dx(2.0_dp switchSpline%c(1) + &
207			3.0_dp * dx * switchSpline%d(1))
208	chrisfen	937	else
209	gezelter	939	call lookupUniformSpline1d(switchSpline, r, sw, dswdr)
210	chrisfen	937	endif
211	gezelter	939
212	gezelter	115	in_switching_region = .true.
213	gezelter	939
214	gezelter	115	return
215			endif
216			else
217			return
218	gezelter	507	endif
219	gezelter	939
220	gezelter	115	end subroutine get_switch
221	chrisfen	725
222	gezelter	115	end module switcheroo