API tools faq
paste
cippo1987

Untitled

cippo1987
Apr 1st, 2019
565
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 9.40 KB | None | 0 0
raw download clone embed print report
  1. subroutine calc_vectors(lattice_11, lattice_12, lattice_13, lattice_21, lattice_22, lattice_23, lattice_31, lattice_32, lattice_33, aa, bb, cc, alpha, beta, gamma, calpha, cbeta,cgamma)  
  2. implicit none
  3. real(kind=8), Intent(In)    :: lattice_11, lattice_12, lattice_13,lattice_21, lattice_22, lattice_23, lattice_31, lattice_32, lattice_33
  4. real(kind=8), Intent(Out)   :: aa, bb, cc
  5. real(kind=8), Intent(Out)   :: alpha, beta, gamma, calpha, cbeta, cgamma
  6. real(kind=8)                                                            :: vol
  7.  
  8.  aa = lattice_11
  9.  bb = sqrt(lattice_21*lattice_21  + lattice_22*lattice_22)
  10.  cc = sqrt(lattice_31*lattice_31 + lattice_32*lattice_32+lattice_33*lattice_33)
  11.  vol = (lattice_31*(lattice_12*lattice_23-lattice_13*lattice_22) + lattice_32*(lattice_13*lattice_21-lattice_11*lattice_23) + lattice_33*(lattice_11*lattice_22-lattice_12*lattice_21))
  12.  
  13. alpha = acos((lattice_21*lattice_31 + lattice_22*lattice_32 + lattice_23*lattice_33)/(bb*cc))
  14. beta  = acos((lattice_11*lattice_31 + lattice_12*lattice_32 + lattice_13*lattice_33)/(aa*cc))
  15. gamma = acos((lattice_11*lattice_21 + lattice_12*lattice_22 + lattice_13*lattice_23)/(bb*aa))
  16.  
  17.  calpha = cos(alpha)
  18.  cbeta  = cos(beta)
  19.  cgamma = cos(gamma)
  20. end subroutine calc_vectors
  21. !!!!
  22. !
  23. !!!!
  24. subroutine read_structure(UnitNum,FileName,num_atom, atom_position)  
  25. implicit none
  26. integer, Intent(In)                                                     :: UnitNum
  27. character (len=*), Intent(In)                                           :: FileName
  28. integer(kind=8)                                                         :: i
  29. integer(kind=8), Intent(In)                                             :: num_atom
  30.  
  31. real(kind=8),    DIMENSION(num_atom, 3), Intent(InOut)                  :: atom_position
  32. !write(*,*)UnitNum, FileName
  33. open(unit=UnitNum, file=FileName, status='old', action='read' )
  34. atom_position = -1
  35. read(1,*)! SKIP   !read position of first specie
  36.  
  37. do i=1, num_atom
  38.     read(1,*)atom_position(i,1), atom_position(i,2), atom_position(i,3)    
  39. end do
  40.  
  41. return
  42. end subroutine read_structure
  43. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  44.  
  45. subroutine rdf(vector_rdf, lung, number_of_couples, couple)
  46. implicit none
  47. integer(kind=8), Intent(InOut)  :: vector_rdf(number_of_couples,600)
  48. real(kind=8), Intent(In)        :: lung
  49. integer(kind=8),Intent(In)      :: number_of_couples, couple
  50. real                    :: bin
  51. integer                 :: position
  52. !write(*,*)lung, "lung"
  53. bin = 0.1
  54. position =  INT(lung/bin) + 1
  55. !write(*,*) "chiamata a rdf: position", position, "bin", bin, "lung", lung
  56. vector_rdf(couple, position)=vector_rdf(couple, position) + 1
  57. RETURN
  58. end subroutine rdf
  59. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  60. subroutine setupvector(vector_rdf,number_of_couples)
  61. implicit none
  62. integer(kind=8), Intent(InOut)      :: vector_rdf(number_of_couples,600)
  63. integer(kind=8),Intent(In)                  :: number_of_couples
  64. vector_rdf=0
  65. RETURN
  66. end subroutine setupvector
  67. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  68. program read_xdatcar
  69. implicit none
  70. integer(kind=8)                                         :: num_structu
  71. integer(kind=8)                                         :: num_atom
  72. integer(kind=8)                                         :: num_species, atoma, atomb
  73. integer(kind=8), DIMENSION(:),    ALLOCATABLE           :: atom_a, atom_b
  74. integer(kind=8)                                         :: number_of_couples
  75.  
  76.   real(kind=8),    DIMENSION(:, :), ALLOCATABLE              :: atom_position
  77. integer(kind=8), DIMENSION(:),    ALLOCATABLE              :: element_quantity
  78. integer(kind=8), DIMENSION(3)                             :: per
  79.   real(kind=8)                                          :: lattice_11, lattice_12, lattice_13
  80. real(kind=8)                                            :: lattice_21, lattice_22, lattice_23
  81. real(kind=8)                                            :: lattice_31, lattice_32, lattice_33
  82.   real(kind=8)                                          :: distance, min_distance, bond
  83.   real(kind=8)                                          :: aa,bb,cc
  84. real(kind=8)                                            :: alpha, beta, gamma
  85. real(kind=8)                                            :: calpha, cbeta, cgamma
  86. integer(kind=8)                                         :: i , ii, jj, kk, jkl, write_i
  87. integer(kind=8)                                         :: unit_output
  88.  
  89. integer(kind=8), DIMENSION(:, :), ALLOCATABLE                           :: rdf_minimo, rdf_all, vector_rdf
  90. !integer(kind=8), DIMENSION(600)                            :: rdf_minimo, rdf_all, vector_rdf
  91.  
  92. real(kind=8)                                            :: x_axis
  93. ! DEFINE FILES
  94. open (unit =  1, file="ONE_XDATCAR", status='old', action='read')  ! trajectory as in HEADER + ALL THE CONFIGURATIONS
  95. open (unit =  2, file="INPUT", status='old', action='read')  ! trajectory as in HEADER + ALL THE CONFIGURATIONS
  96.  
  97. per(1) = -1
  98. per(2) = 0
  99. per(3) = 1
  100. read(1,*) num_structu, num_species
  101. read(1,*) !skip
  102. read(1,*) lattice_11, lattice_12, lattice_13
  103. read(1,*) lattice_21, lattice_22, lattice_23
  104. read(1,*) lattice_31, lattice_32, lattice_33
  105. read(1,*) ! skip elements name  
  106.  
  107. !READ INPUT
  108. read(2,*)number_of_couples
  109.  
  110. ALLOCATE(rdf_minimo(number_of_couples,600))
  111. ALLOCATE(rdf_all(number_of_couples,600))
  112. ALLOCATE(vector_rdf(number_of_couples,600))
  113. ALLOCATE(atom_a(number_of_couples))
  114. ALLOCATE(atom_b(number_of_couples))
  115. atom_a = 0
  116. atom_b = 0
  117.  
  118. do jkl=1, number_of_couples
  119. read(2,*)atom_a(jkl),atom_b(jkl)
  120. end do
  121.  
  122.  
  123.  
  124. rdf_minimo = 0
  125. rdf_all = 0
  126. !RESET VECTOR RDF
  127. call setupvector(vector_rdf,number_of_couples)
  128. call calc_vectors(lattice_11, lattice_12, lattice_13, lattice_21, lattice_22, lattice_23, lattice_31, lattice_32, lattice_33, aa, bb, cc, alpha, beta, gamma, calpha, cbeta,cgamma)
  129. ! THIS SUB calculates the lenght of the cell and the angles.
  130. !write(*,*)calpha
  131. ALLOCATE(element_quantity(num_species))
  132. element_quantity = 0
  133. SELECT CASE (num_species) ! I don't know how many elements I do have.
  134.    CASE (1)
  135. read(1,*) element_quantity(1)
  136.    CASE (2)
  137. read(1,*) element_quantity(1), element_quantity(2)
  138.    CASE (3)
  139. read(1,*) element_quantity(1), element_quantity(2), element_quantity(3)
  140.    CASE (4)
  141. read(1,*) element_quantity(1), element_quantity(2), element_quantity(3), element_quantity(4)
  142.    CASE (5)
  143. read(1,*) element_quantity(1), element_quantity(2), element_quantity(3), element_quantity(4), element_quantity(5)
  144.    CASE (6)
  145. read(1,*) element_quantity(1), element_quantity(2), element_quantity(3), element_quantity(4), element_quantity(5), element_quantity(6)
  146.    CASE (7)
  147. read(1,*) element_quantity(1), element_quantity(2), element_quantity(3), element_quantity(4), element_quantity(5), element_quantity(6), element_quantity(7)
  148.    CASE (8)
  149. read(1,*) element_quantity(1), element_quantity(2), element_quantity(3), element_quantity(4), element_quantity(5), element_quantity(6), element_quantity(7), element_quantity(8)
  150. END SELECT
  151. !write(*,*)'num species', num_species
  152. num_atom = 0 ! RESET variable
  153. do i=1,num_species
  154.     num_atom = num_atom + element_quantity(i)
  155. end do
  156. !write(*,*)"num atom", num_atom
  157. ALLOCATE(atom_position(num_atom,3) )
  158. atom_position = 0
  159. !!!!!!!!!!!!!!!!!!
  160. !!!!!!!!!!!!!!!!!!
  161. !!!    START   !!!
  162. !!!!!!!!!!!!!!!!!!
  163. !!!!!!!!!!!!!!!!!!
  164. do i=1,num_structu
  165.  
  166.     call read_structure(1,"ONE_XDATCAR",num_atom, atom_position) ! 1 is the file unit... !ONE_XDATCAR is the file name
  167.     !write(*,*)distance(1,2, aa, bb, cc, calpha, cbeta, cgamma, atom_position, num_atom)
  168.     ! After one structure is read, calculate d, print it.
  169.     do jkl=1,number_of_couples
  170.         atoma = atom_a(jkl)
  171.                 atomb = atom_b(jkl)
  172.  
  173.         min_distance = 10000.45   !Necessary to void random number
  174.             do ii=1,3
  175.                 do jj=1,3
  176.                     do kk=1,3
  177.                        
  178.                         atom_position(atomb,1) = (atom_position(atomb,1) + per(ii))     ! generate all the 27 possible distance. Not very clever.
  179.                         atom_position(atomb,2) = (atom_position(atomb,2) + per(jj))
  180.                         atom_position(atomb,3) = (atom_position(atomb,3) + per(kk))
  181.                        
  182.                         bond=distance(atoma, atomb, aa, bb, cc, calpha, cbeta, cgamma, atom_position, num_atom)
  183.                         if ( (bond.lt.min_distance) ) then      !.and.(pbclung.gt.0)) then  !write(*,*)"MIN", bond, min_distance
  184.                             min_distance = bond
  185.                         end if  !write(*,*)"this is bond, which is passed as lung", bond
  186.                         !subroutine rdf(vector_rdf, lung, number_of_couples, couple)
  187.                         call rdf(rdf_all, bond, number_of_couples, jkl)                 !write(*,*)"b", bond
  188.                        
  189.                         atom_position(atomb,1) = (atom_position(atomb,1) - per(ii))     ! generate all the 27 possible distance. Not very clever.
  190.                         atom_position(atomb,2) = (atom_position(atomb,2) - per(jj))
  191.                         atom_position(atomb,3) = (atom_position(atomb,3) - per(kk))
  192.  
  193.  
  194.                     end do  !ii
  195.                 end do  !! jj
  196.             end do !! kk                !write(*,*)"minima",min_distance
  197.         call rdf(rdf_minimo, min_distance, number_of_couples, jkl)                  !write(*,*)"this is min_distance, which is passed as lung", min_distance
  198.     end do   ! jkl
  199.    
  200. end do ! num_structure
  201.  
  202.     do jkl=1, number_of_couples
  203.         unit_output = 10 + jkl
  204.         x_axis = 0.0
  205.         do write_i=1,600
  206.             write(unit_output,*)x_axis, rdf_minimo(jkl,write_i)
  207.             write(unit_output+30,*)x_axis, rdf_all(jkl,write_i)
  208.             x_axis = x_axis + 0.1
  209.         end do ! jkl writing files
  210.     end do ! i writing files
  211.  
  212. end program read_xdatcar
  213. !!!!!!!!!BEGIN FUNCTION
  214.  
  215. real(kind=8) function distance(atom1, atom2, aa, bb, cc, calpha, cbeta, cgamma, atom_position, num_atom)
  216. implicit none
  217.  
  218. integer(kind=8),      intent(in)                                                :: atom1, atom2, num_atom
  219. real(kind=8), intent(in)                                                        :: aa, bb, cc, calpha, cbeta, cgamma
  220. real(kind=8),    DIMENSION(num_atom,3), Intent(In)                              :: atom_position
  221. real(kind=8)                                                                    :: dfx, dfy, dfz
  222.  dfx = (atom_position(atom1,1) - atom_position(atom2,1))
  223. dfy = (atom_position(atom1,2) - atom_position(atom2,2))
  224. dfz = (atom_position(atom1,3) - atom_position(atom2,3))
  225.  distance = sqrt(aa*aa*dfx*dfx + bb*bb*dfy*dfy + cc*cc*dfz*dfz + 2*bb*cc*calpha*dfy*dfz + 2*cc*aa*cbeta*dfz*dfx + 2*aa*bb*cgamma*dfx*dfy)
  226.  RETURN  !no need of return if contains is used
  227. end function distance
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!