From e76187bd3c177ce8ecc536efc68210fa852f137d Mon Sep 17 00:00:00 2001
From: Alex <aselimov3@gatech.edu>
Date: Mon, 27 Jan 2020 10:16:07 -0500
Subject: [PATCH] Square loop option

---
 src/opt_disl.f90 | 75 +++++++++++++++++++++++++++++-------------------
 1 file changed, 45 insertions(+), 30 deletions(-)

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