Statistics
| Revision:

## root / src / Calc_Xprim.f90 @ 7

 1 SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)  !  ! This subroutine uses the description of a list of Coordinates  ! to compute the values of the coordinates for a given geometry.  !  !!!!!!!!!!  ! Input:  ! Na: INTEGER, Number of atoms  ! x,y,z(Na): REAL, cartesian coordinates of the considered geometry  ! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes  ! NPrim, INTEGER: Number of coordinates to compute  !  ! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for  ! a former geometry. Useful for Dihedral angles evolution...  !!!!!!!!!!!  ! Output:  ! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates  !  !!!!!!!!!   Use VarTypes   Use Io_module   Use Path_module, only : pi   IMPLICIT NONE   Type (ListCoord), POINTER :: Coordinate   INTEGER(KINT), INTENT(IN) :: Nat,NPrim   REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)   REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim)   REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)   Type (ListCoord), POINTER :: ScanCoord   real(KREAL) :: vx1,vy1,vz1,norm1   real(KREAL) :: vx2,vy2,vz2,norm2   real(KREAL) :: vx3,vy3,vz3,norm3   real(KREAL) :: vx4,vy4,vz4,norm4   real(KREAL) :: vx5,vy5,vz5,norm5   INTEGER(KINT) :: I   REAL(KREAL) :: DiheTmp   LOGICAL :: debug, debugPFL   INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID   FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)   use Path_module, only : Pi,KINT, KREAL   real(KREAL) :: v1x,v1y,v1z,norm1   real(KREAL) :: v2x,v2y,v2z,norm2   real(KREAL) :: angle   END FUNCTION ANGLE   FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)   use Path_module, only : Pi,KINT, KREAL   real(KREAL) :: v1x,v1y,v1z,norm1   real(KREAL) :: v2x,v2y,v2z,norm2   real(KREAL) :: SinAngle   END FUNCTION SINANGLE   FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)   use Path_module, only : Pi,KINT, KREAL   real(KREAL) :: v1x,v1y,v1z,norm1   real(KREAL) :: v2x,v2y,v2z,norm2   real(KREAL) :: v3x,v3y,v3z,norm3   real(KREAL) :: angle_d,ca,sa   END FUNCTION ANGLE_D   END INTERFACE   debug=valid("Calc_Xprim")   debugPFL=valid("BakerPFL")   if (debug) Call Header("Entering Cal_XPrim")   IF (debug) THEN   WRITE(*,*) "DBG Calc_Xprim, geom"   DO I=1,Nat   WRITE(*,'(1X,I5,3(1X,F15.3))') I,X(I),y(i),z(i)   END DO   IF (Present(XPrimRef)) THEN   WRITE(*,*) "XPrimRef"   WRITE(*,'(15(1X,F10.6))') XPrimRef   END IF   WRITE(*,*) "NPrim:",NPrim   END IF   WRITE(*,*) "Coordinate:",associated(Coordinate),coordinate%Type   ScanCoord => Coordinate   WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type  ! WRITE(*,*) "coucou"   I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).   DO WHILE (Associated(ScanCoord%next))   I=I+1  ! WRITE(*,*) i   SELECT CASE (ScanCoord%Type)   CASE ('BOND')   Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)   Xprimitive(I) = Norm2   CASE ('ANGLE')   Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1)   Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)   Xprimitive(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.   CASE ('DIHEDRAL')   Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)   Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)   Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)   Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &   vx4,vy4,vz4,norm4)   Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &   vx5,vy5,vz5,norm5)     DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &   vx2,vy2,vz2,norm2)   Xprimitive(I) = DiheTmp*Pi/180.  ! PFL 15th March 2008 ->  ! I think that the test on changes less than Pi should be enough  !! We treat large dihedral angles differently as +180=-180 mathematically and physically  !! but this causes lots of troubles when converting baker to cart.  !! So we ensure that large dihedral angles always have the same sign  ! if (abs(ScanCoord%SignDihedral).NE.1) THEN  ! ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp))  ! ELSE  ! If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN  ! Xprimitive(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi  ! END IF  ! END IF  !!!! <- PFL 15th March 2008  ! Another test... will all this be consistent ???  ! We want the shortest path, so we check that the change in dihedral angles is less than Pi:   IF (Present(XPrimRef)) THEN   IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN   if (debug) THEN   WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I)   WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi   END IF   if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN   Xprimitive(I)=Xprimitive(I)-2.d0*Pi   ELSE   Xprimitive(I)=Xprimitive(I)+2.d0*Pi   END IF   END IF   if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi   END IF   CASE DEFAULT   WRITE(*,*) "Scancoord%type unknown for I",I,scancoord%type   END SELECT   ScanCoord => ScanCoord%next   END DO ! matches DO WHILE (Associated(ScanCoord%next))   IF (debug) THEN   WRITE(*,*) "DBG Calc_Xprim Values"   WRITE(*,*) "Found ",I," primitives"     ScanCoord=>Coordinate   I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).   DO WHILE (Associated(ScanCoord%next))   I=I+1   SELECT CASE (ScanCoord%Type)   CASE ('BOND')   WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(I)   CASE ('ANGLE')   WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &   ScanCoord%At2, ScanCoord%At3,Xprimitive(I)*180./Pi   CASE ('DIHEDRAL')   WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&   ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(I)*180./Pi   END SELECT   ScanCoord => ScanCoord%next   END DO ! matches DO WHILE (Associated(ScanCoord%next))   END IF   if (debug) Call Header(" Cal_XPrim OVER ")  END SUBROUTINE Calc_Xprim