diff --git a/src/opt_disl.f90 b/src/opt_disl.f90 index 4510ede..00886e5 100644 --- a/src/opt_disl.f90 +++ b/src/opt_disl.f90 @@ -291,36 +291,51 @@ module opt_disl a3 = 3 end select - !Calculate loop perimeter - perimeter = 2.0_dp*pi*loop_radius - - !Define the number of points forming the loop - ! The following criteria are used as a trade-off between - ! good accuracy and computational efficiency: - ! - each dislocation segment should have a length of 5 angströms; - ! - the loop should contain at least 3 points - ! (for very small loops, this will result in segments shorter than 5 A); - ! - there should not be more than 100 points - ! (for very large loops, this will result in segments longer than 5 A). - Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) ) - - !angle between two consecutive points - theta = 2.0_dp*pi / dble(Npoints) - - !allocate xLoop - allocate(xLoop(Npoints,3)) - xLoop(:,:) = 0.0_dp - - !Calculate the position of each point in the loop - angle = 0.0_dp - do i = 1, size(xLoop,1) - xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle) - xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle) - xLoop(i,a3) = centroid(a3) - ! Increment angle for next point - angle = angle + theta - end do - + if(loop_radius < 0.0_dp) then + xLoop(:,:) = 0.d0 + xLoop(1,a1) = centroid(1) - loop_radius + xLoop(1,a2) = centroid(2) - loop_radius + xLoop(1,a3) = centroid(3) + xLoop(2,a1) = centroid(1) + loop_radius + xLoop(2,a2) = centroid(2) - loop_radius + xLoop(2,a3) = centroid(3) + xLoop(3,a1) = centroid(1) + loop_radius + xLoop(3,a2) = centroid(2) + loop_radius + xLoop(3,a3) = centroid(3) + xLoop(4,a1) = centroid(1) - loop_radius + xLoop(4,a2) = centroid(2) + loop_radius + xLoop(4,a3) = centroid(3) + else + !Calculate loop perimeter + perimeter = 2.0_dp*pi*loop_radius + + !Define the number of points forming the loop + ! The following criteria are used as a trade-off between + ! good accuracy and computational efficiency: + ! - each dislocation segment should have a length of 5 angströms; + ! - the loop should contain at least 3 points + ! (for very small loops, this will result in segments shorter than 5 A); + ! - there should not be more than 100 points + ! (for very large loops, this will result in segments longer than 5 A). + Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) ) + + !angle between two consecutive points + theta = 2.0_dp*pi / dble(Npoints) + + !allocate xLoop + allocate(xLoop(Npoints,3)) + xLoop(:,:) = 0.0_dp + + !Calculate the position of each point in the loop + angle = 0.0_dp + do i = 1, size(xLoop,1) + xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle) + xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle) + xLoop(i,a3) = centroid(a3) + ! Increment angle for next point + angle = angle + theta + end do + end if !Now actually calculate the displacement created by a loop for every atom do i = 1, atom_num u = 0.0_dp