ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/suttonchen.F90
Revision: 2427
Committed: Mon Nov 14 21:29:36 2005 UTC (18 years, 9 months ago) by chuckv
File size: 19084 byte(s)
Log Message:
Sutton-Chen almost done. Just need to fix do_forces to use Sutton-Chen.

File Contents

# Content
1 !!
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 !! 1. Acknowledgement of the program authors must be made in any
10 !! publication of scientific results based in part on use of the
11 !! program. An acceptable form of acknowledgement is citation of
12 !! the article in which the program was described (Matthew
13 !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 !! Parallel Simulation Engine for Molecular Dynamics,"
16 !! J. Comput. Chem. 26, pp. 252-271 (2005))
17 !!
18 !! 2. Redistributions of source code must retain the above copyright
19 !! notice, this list of conditions and the following disclaimer.
20 !!
21 !! 3. Redistributions in binary form must reproduce the above copyright
22 !! notice, this list of conditions and the following disclaimer in the
23 !! documentation and/or other materials provided with the
24 !! distribution.
25 !!
26 !! This software is provided "AS IS," without a warranty of any
27 !! kind. All express or implied conditions, representations and
28 !! warranties, including any implied warranty of merchantability,
29 !! fitness for a particular purpose or non-infringement, are hereby
30 !! excluded. The University of Notre Dame and its licensors shall not
31 !! be liable for any damages suffered by licensee as a result of
32 !! using, modifying or distributing the software or its
33 !! derivatives. In no event will the University of Notre Dame or its
34 !! licensors be liable for any lost revenue, profit or data, or for
35 !! direct, indirect, special, consequential, incidental or punitive
36 !! damages, however caused and regardless of the theory of liability,
37 !! arising out of the use of or inability to use software, even if the
38 !! University of Notre Dame has been advised of the possibility of
39 !! such damages.
40 !!
41
42 !! Impliments Sutton-Chen Metallic Potential
43 !! See A.P.SUTTON and J.CHEN,PHIL MAG LETT 61,139-146,1990
44
45
46 module suttonchen
47 use simulation
48 use force_globals
49 use status
50 use atype_module
51 use vector_class
52 #ifdef IS_MPI
53 use mpiSimulation
54 #endif
55 implicit none
56 PRIVATE
57 #define __FORTRAN90
58 #include "UseTheForce/DarkSide/fInteractionMap.h"
59
60 INTEGER, PARAMETER :: DP = selected_real_kind(15)
61
62 logical, save :: SC_FF_initialized = .false.
63 integer, save :: SC_Mixing_Policy
64 real(kind = dp), save :: SC_rcut
65 logical, save :: haveRcut = .false.
66 logical, save :: haveMixingMap = .false.
67 logical, save :: useGeometricDistanceMixing = .false.
68
69
70
71
72 character(len = statusMsgSize) :: errMesg
73 integer :: sc_err
74
75 character(len = 200) :: errMsg
76 character(len=*), parameter :: RoutineName = "Sutton-Chen MODULE"
77 !! Logical that determines if eam arrays should be zeroed
78 logical :: cleanme = .true.
79 logical :: nmflag = .false.
80
81
82 type, private :: SCtype
83 integer :: atid
84 real(kind=dp) :: c
85 real(kind=dp) :: m
86 real(kind=dp) :: n
87 real(kind=dp) :: alpha
88 real(kind=dp) :: epsilon
89 real(kind=dp) :: sc_rcut
90 end type SCtype
91
92
93 !! Arrays for derivatives used in force calculation
94 real( kind = dp), dimension(:), allocatable :: rho
95 real( kind = dp), dimension(:), allocatable :: frho
96 real( kind = dp), dimension(:), allocatable :: dfrhodrho
97
98
99
100 !! Arrays for MPI storage
101 #ifdef IS_MPI
102 real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
103 real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
104 real( kind = dp),save, dimension(:), allocatable :: frho_row
105 real( kind = dp),save, dimension(:), allocatable :: frho_col
106 real( kind = dp),save, dimension(:), allocatable :: rho_row
107 real( kind = dp),save, dimension(:), allocatable :: rho_col
108 real( kind = dp),save, dimension(:), allocatable :: rho_tmp
109 #endif
110
111 type, private :: SCTypeList
112 integer :: nSCTypes = 0
113 integer :: currentSCtype = 0
114
115 type (SCtype), pointer :: SCtypes(:) => null()
116 integer, pointer :: atidToSCtype(:) => null()
117 end type SCTypeList
118
119 type (SCTypeList), save :: SCList
120
121
122
123
124 type :: MixParameters
125 real(kind=DP) :: alpha
126 real(kind=DP) :: epsilon
127 real(kind=DP) :: m
128 real(Kind=DP) :: n
129 real(kind=DP) :: vpair_pot
130 real(kind=dp) :: rCut
131 logical :: rCutWasSet = .false.
132
133 end type MixParameters
134
135 type(MixParameters), dimension(:,:), allocatable :: MixingMap
136
137
138
139 public :: setCutoffSC
140 public :: do_SC_pair
141 public :: newSCtype
142 public :: calc_SC_prepair_rho
143 public :: clean_SC
144 public :: destroySCtypes
145 public :: getSCCut
146 ! public :: setSCDefaultCutoff
147 ! public :: setSCUniformCutoff
148 public :: useGeometricMixing
149
150 contains
151
152
153 subroutine newSCtype(c_ident,c,m,n,alpha,epsilon,status)
154 real (kind = dp ) :: c ! Density Scaling
155 real (kind = dp ) :: m ! Density Exponent
156 real (kind = dp ) :: n ! Pair Potential Exponent
157 real (kind = dp ) :: alpha ! Length Scaling
158 real (kind = dp ) :: epsilon ! Energy Scaling
159
160
161 integer :: c_ident
162 integer :: status
163
164 integer :: nAtypes,nSCTypes,myATID
165 integer :: maxVals
166 integer :: alloc_stat
167 integer :: current
168 integer,pointer :: Matchlist(:) => null()
169
170 status = 0
171
172
173 !! Assume that atypes has already been set and get the total number of types in atypes
174
175
176 ! check to see if this is the first time into
177 if (.not.associated(SCList%SCTypes)) then
178 call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
179 SCList%nSCtypes = nSCtypes
180 allocate(SCList%SCTypes(nSCTypes))
181 nAtypes = getSize(atypes)
182 allocate(SCList%atidToSCType(nAtypes))
183 end if
184
185 SCList%currentSCType = SCList%currentSCType + 1
186 current = SCList%currentSCType
187
188 myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
189 SCList%atidToSCType(myATID) = current
190
191
192
193 SCList%SCTypes(current)%atid = c_ident
194 SCList%SCTypes(current)%alpha = alpha
195 SCList%SCTypes(current)%c = c
196 SCList%SCTypes(current)%m = m
197 SCList%SCTypes(current)%n = n
198 SCList%SCTypes(current)%epsilon = epsilon
199 end subroutine newSCtype
200
201
202 subroutine destroySCTypes()
203 if (associated(SCList%SCtypes)) then
204 deallocate(SCList%SCTypes)
205 SCList%SCTypes=>null()
206 end if
207 if (associated(SCList%atidToSCtype)) then
208 deallocate(SCList%atidToSCtype)
209 SCList%atidToSCtype=>null()
210 end if
211
212
213 end subroutine destroySCTypes
214
215
216
217 function getSCCut(atomID) result(cutValue)
218 integer, intent(in) :: atomID
219 integer :: scID
220 real(kind=dp) :: cutValue
221
222 scID = SCList%atidToSCType(atomID)
223 cutValue = SCList%SCTypes(scID)%sc_rcut
224 end function getSCCut
225
226
227
228
229 subroutine createMixingMap()
230 integer :: nSCtypes, i, j
231 real ( kind = dp ) :: e1, e2,m1,m2,alpha1,alpha2,n1,n2
232 real ( kind = dp ) :: rcut6, tp6, tp12
233 logical :: isSoftCore1, isSoftCore2, doShift
234
235 if (SCList%currentSCtype == 0) then
236 call handleError("SuttonChen", "No members in SCMap")
237 return
238 end if
239
240 nSCtypes = SCList%nSCtypes
241
242 if (.not. allocated(MixingMap)) then
243 allocate(MixingMap(nSCtypes, nSCtypes))
244 endif
245
246 do i = 1, nSCtypes
247
248 e1 = SCList%SCtypes(i)%epsilon
249 m1 = SCList%SCtypes(i)%m
250 n1 = SCList%SCtypes(i)%n
251 alpha1 = SCList%SCtypes(i)%alpha
252
253 do j = i, nSCtypes
254
255
256 e2 = SCList%SCtypes(j)%epsilon
257 m2 = SCList%SCtypes(j)%m
258 n2 = SCList%SCtypes(j)%n
259 alpha2 = SCList%SCtypes(j)%alpha
260
261 if (useGeometricDistanceMixing) then
262 MixingMap(i,j)%alpha = sqrt(alpha1 * alpha2) !SC formulation
263 else
264 MixingMap(i,j)%alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
265 endif
266
267 MixingMap(i,j)%epsilon = sqrt(e1 * e2)
268 MixingMap(i,j)%m = 0.5_dp*(m1+m2)
269 MixingMap(i,j)%n = 0.5_dp*(n1+n2)
270 MixingMap(i,j)%alpha = 0.5_dp*(alpha1+alpha2)
271 MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
272 MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
273 (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
274 if (i.ne.j) then
275 MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
276 MixingMap(j,i)%m = MixingMap(i,j)%m
277 MixingMap(j,i)%n = MixingMap(i,j)%n
278 MixingMap(j,i)%alpha = MixingMap(i,j)%alpha
279 MixingMap(j,i)%rcut = MixingMap(i,j)%rcut
280 MixingMap(j,i)%vpair_pot = MixingMap(i,j)%vpair_pot
281 endif
282 enddo
283 enddo
284
285 haveMixingMap = .true.
286
287 end subroutine createMixingMap
288
289
290 !! routine checks to see if array is allocated, deallocates array if allocated
291 !! and then creates the array to the required size
292 subroutine allocateSC(status)
293 integer, intent(out) :: status
294
295 #ifdef IS_MPI
296 integer :: nAtomsInRow
297 integer :: nAtomsInCol
298 #endif
299 integer :: alloc_stat
300
301
302 status = 0
303 #ifdef IS_MPI
304 nAtomsInRow = getNatomsInRow(plan_atom_row)
305 nAtomsInCol = getNatomsInCol(plan_atom_col)
306 #endif
307
308
309
310 if (allocated(frho)) deallocate(frho)
311 allocate(frho(nlocal),stat=alloc_stat)
312 if (alloc_stat /= 0) then
313 status = -1
314 return
315 end if
316
317 if (allocated(rho)) deallocate(rho)
318 allocate(rho(nlocal),stat=alloc_stat)
319 if (alloc_stat /= 0) then
320 status = -1
321 return
322 end if
323
324 if (allocated(dfrhodrho)) deallocate(dfrhodrho)
325 allocate(dfrhodrho(nlocal),stat=alloc_stat)
326 if (alloc_stat /= 0) then
327 status = -1
328 return
329 end if
330
331 #ifdef IS_MPI
332
333 if (allocated(rho_tmp)) deallocate(rho_tmp)
334 allocate(rho_tmp(nlocal),stat=alloc_stat)
335 if (alloc_stat /= 0) then
336 status = -1
337 return
338 end if
339
340
341 if (allocated(frho_row)) deallocate(frho_row)
342 allocate(frho_row(nAtomsInRow),stat=alloc_stat)
343 if (alloc_stat /= 0) then
344 status = -1
345 return
346 end if
347 if (allocated(rho_row)) deallocate(rho_row)
348 allocate(rho_row(nAtomsInRow),stat=alloc_stat)
349 if (alloc_stat /= 0) then
350 status = -1
351 return
352 end if
353 if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
354 allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
355 if (alloc_stat /= 0) then
356 status = -1
357 return
358 end if
359
360
361 ! Now do column arrays
362
363 if (allocated(frho_col)) deallocate(frho_col)
364 allocate(frho_col(nAtomsInCol),stat=alloc_stat)
365 if (alloc_stat /= 0) then
366 status = -1
367 return
368 end if
369 if (allocated(rho_col)) deallocate(rho_col)
370 allocate(rho_col(nAtomsInCol),stat=alloc_stat)
371 if (alloc_stat /= 0) then
372 status = -1
373 return
374 end if
375 if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
376 allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat)
377 if (alloc_stat /= 0) then
378 status = -1
379 return
380 end if
381
382 #endif
383
384 end subroutine allocateSC
385
386 !! C sets rcut to be the largest cutoff of any atype
387 !! present in this simulation. Doesn't include all atypes
388 !! sim knows about, just those in the simulation.
389 subroutine setCutoffSC(rcut, status)
390 real(kind=dp) :: rcut
391 integer :: status
392 status = 0
393
394 SC_rcut = rcut
395
396 end subroutine setCutoffSC
397
398 subroutine useGeometricMixing()
399 useGeometricDistanceMixing = .true.
400 haveMixingMap = .false.
401 return
402 end subroutine useGeometricMixing
403
404
405
406
407
408
409
410
411
412 subroutine clean_SC()
413
414 ! clean non-IS_MPI first
415 frho = 0.0_dp
416 rho = 0.0_dp
417 dfrhodrho = 0.0_dp
418 ! clean MPI if needed
419 #ifdef IS_MPI
420 frho_row = 0.0_dp
421 frho_col = 0.0_dp
422 rho_row = 0.0_dp
423 rho_col = 0.0_dp
424 rho_tmp = 0.0_dp
425 dfrhodrho_row = 0.0_dp
426 dfrhodrho_col = 0.0_dp
427 #endif
428 end subroutine clean_SC
429
430
431
432 !! Calculates rho_r
433 subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq)
434 integer :: atom1,atom2
435 real(kind = dp), dimension(3) :: d
436 real(kind = dp), intent(inout) :: r
437 real(kind = dp), intent(inout) :: rijsq
438 ! value of electron density rho do to atom i at atom j
439 real(kind = dp) :: rho_i_at_j
440 ! value of electron density rho do to atom j at atom i
441 real(kind = dp) :: rho_j_at_i
442
443 ! we don't use the derivatives, dummy variables
444 real( kind = dp) :: drho
445 integer :: sc_err
446
447 integer :: atid1,atid2 ! Global atid
448 integer :: myid_atom1 ! EAM atid
449 integer :: myid_atom2
450
451
452 ! check to see if we need to be cleaned at the start of a force loop
453
454
455
456
457 #ifdef IS_MPI
458 Atid1 = Atid_row(Atom1)
459 Atid2 = Atid_col(Atom2)
460 #else
461 Atid1 = Atid(Atom1)
462 Atid2 = Atid(Atom2)
463 #endif
464
465 Myid_atom1 = SCList%atidtoSCtype(Atid1)
466 Myid_atom2 = SCList%atidtoSCtype(Atid2)
467
468
469
470 rho_i_at_j = (MixingMap(Myid_atom1,Myid_atom2)%alpha/r)&
471 **MixingMap(Myid_atom1,Myid_atom2)%m
472 rho_j_at_i = rho_i_at_j
473
474 #ifdef IS_MPI
475 rho_col(atom2) = rho_col(atom2) + rho_i_at_j
476 rho_row(atom1) = rho_row(atom1) + rho_j_at_i
477 #else
478 rho(atom2) = rho(atom2) + rho_i_at_j
479 rho(atom1) = rho(atom1) + rho_j_at_i
480 #endif
481
482 end subroutine calc_sc_prepair_rho
483
484
485 !! Calculate the rho_a for all local atoms
486 subroutine calc_sc_preforce_Frho(nlocal,pot)
487 integer :: nlocal
488 real(kind=dp) :: pot
489 integer :: i,j
490 integer :: atom
491 real(kind=dp) :: U,U1,U2
492 integer :: atype1
493 integer :: atid1
494 integer :: myid
495
496
497 cleanme = .true.
498 !! Scatter the electron density from pre-pair calculation back to local atoms
499 #ifdef IS_MPI
500 call scatter(rho_row,rho,plan_atom_row,sc_err)
501 if (sc_err /= 0 ) then
502 write(errMsg,*) " Error scattering rho_row into rho"
503 call handleError(RoutineName,errMesg)
504 endif
505 call scatter(rho_col,rho_tmp,plan_atom_col,sc_err)
506 if (sc_err /= 0 ) then
507 write(errMsg,*) " Error scattering rho_col into rho"
508 call handleError(RoutineName,errMesg)
509
510 endif
511
512 rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
513 #endif
514
515
516
517 !! Calculate F(rho) and derivative
518 do atom = 1, nlocal
519 Myid = SCList%atidtoSctype(Atid(atom))
520 frho(atom) = -SCList%SCTypes(Myid)%c * &
521 SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
522
523 dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
524 pot = pot + u
525 enddo
526
527 #ifdef IS_MPI
528 !! communicate f(rho) and derivatives back into row and column arrays
529 call gather(frho,frho_row,plan_atom_row, sc_err)
530 if (sc_err /= 0) then
531 call handleError("cal_eam_forces()","MPI gather frho_row failure")
532 endif
533 call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, sc_err)
534 if (sc_err /= 0) then
535 call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
536 endif
537 call gather(frho,frho_col,plan_atom_col, sc_err)
538 if (sc_err /= 0) then
539 call handleError("cal_eam_forces()","MPI gather frho_col failure")
540 endif
541 call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, sc_err)
542 if (sc_err /= 0) then
543 call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
544 endif
545 #endif
546
547
548 end subroutine calc_sc_preforce_Frho
549
550
551 !! Does Sutton-Chen pairwise Force calculation.
552 subroutine do_sc_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
553 pot, f, do_pot)
554 !Arguments
555 integer, intent(in) :: atom1, atom2
556 real( kind = dp ), intent(in) :: rij, r2
557 real( kind = dp ) :: pot, sw, vpair
558 real( kind = dp ), dimension(3,nLocal) :: f
559 real( kind = dp ), intent(in), dimension(3) :: d
560 real( kind = dp ), intent(inout), dimension(3) :: fpair
561
562 logical, intent(in) :: do_pot
563
564 real( kind = dp ) :: drdx,drdy,drdz
565 real( kind = dp ) :: dvpdr
566 real( kind = dp ) :: drhodr
567 real( kind = dp ) :: dudr
568 real( kind = dp ) :: rcij
569 real( kind = dp ) :: drhoidr,drhojdr
570 real( kind = dp ) :: Fx,Fy,Fz
571 real( kind = dp ) :: r,d2pha,phb,d2phb
572 real( kind = dp ) :: pot_temp,vptmp
573 real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
574 integer :: id1,id2
575 integer :: mytype_atom1
576 integer :: mytype_atom2
577 integer :: atid1,atid2
578 !Local Variables
579
580 ! write(*,*) "Frho: ", Frho(atom1)
581
582
583 dvpdr = 0.0E0_DP
584
585
586 #ifdef IS_MPI
587 atid1 = atid_row(atom1)
588 atid2 = atid_col(atom2)
589 #else
590 atid1 = atid(atom1)
591 atid2 = atid(atom2)
592 #endif
593
594 mytype_atom1 = SCList%atidToSCType(atid1)
595 mytype_atom2 = SCList%atidTOSCType(atid2)
596
597 drdx = d(1)/rij
598 drdy = d(2)/rij
599 drdz = d(3)/rij
600
601
602 epsilonij = MixingMap(mytype_atom1,mytype_atom2)%epsilon
603 aij = MixingMap(mytype_atom1,mytype_atom2)%alpha
604 nij = MixingMap(mytype_atom1,mytype_atom2)%n
605 mij = MixingMap(mytype_atom1,mytype_atom2)%m
606 vcij = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot
607
608 vptmp = epsilonij*((aij/r)**nij)
609
610
611 dvpdr = -nij*vptmp/r
612 drhodr = -mij*((aij/r)**mij)/r
613
614
615 dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
616 + dvpdr
617
618 pot_temp = vptmp + vcij
619
620
621 #ifdef IS_MPI
622 dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) &
623 + dvpdr
624
625 #else
626 dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
627 + dvpdr
628 #endif
629
630
631 fx = dudr * drdx
632 fy = dudr * drdy
633 fz = dudr * drdz
634
635
636 #ifdef IS_MPI
637 if (do_pot) then
638 pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + (pot_temp)*0.5
639 pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + (pot_temp)*0.5
640 end if
641
642 f_Row(1,atom1) = f_Row(1,atom1) + fx
643 f_Row(2,atom1) = f_Row(2,atom1) + fy
644 f_Row(3,atom1) = f_Row(3,atom1) + fz
645
646 f_Col(1,atom2) = f_Col(1,atom2) - fx
647 f_Col(2,atom2) = f_Col(2,atom2) - fy
648 f_Col(3,atom2) = f_Col(3,atom2) - fz
649 #else
650
651 if(do_pot) then
652 pot = pot + pot_temp
653 end if
654
655 f(1,atom1) = f(1,atom1) + fx
656 f(2,atom1) = f(2,atom1) + fy
657 f(3,atom1) = f(3,atom1) + fz
658
659 f(1,atom2) = f(1,atom2) - fx
660 f(2,atom2) = f(2,atom2) - fy
661 f(3,atom2) = f(3,atom2) - fz
662 #endif
663
664
665 #ifdef IS_MPI
666 id1 = AtomRowToGlobal(atom1)
667 id2 = AtomColToGlobal(atom2)
668 #else
669 id1 = atom1
670 id2 = atom2
671 #endif
672
673 if (molMembershipList(id1) .ne. molMembershipList(id2)) then
674
675 fpair(1) = fpair(1) + fx
676 fpair(2) = fpair(2) + fy
677 fpair(3) = fpair(3) + fz
678
679 endif
680
681
682 end subroutine do_sc_pair
683
684
685
686 end module suttonchen