Commit 6adc4523 authored by ofgeorg's avatar ofgeorg

initial submit to the new SVN

parent 69875d4c
SUBROUTINE LJ_POT_COMP(res,vars,m,n,lbounds,ubounds)
USE cmaes_param_mod
USE cmaesOpts_mod
!Parameters
REAL(MK), DIMENSION(n), INTENT(out) :: res
REAL(MK), DIMENSION(m,n),INTENT(in) :: vars
INTEGER, INTENT(in) :: m ! (3*LJ_N-6) dim vector
INTEGER, INTENT(in) :: n
REAL(MK), OPTIONAL :: lbounds
REAL(MK), OPTIONAL :: ubounds
! Locally used variables
REAL(MK) :: dist
REAL(MK), DIMENSION((m+6)/3,3) :: part_mat
REAL(MK), DIMENSION(1,3) :: part_mean
REAL(MK) :: ub,lb
REAL(MK) :: LJsigma,LJeps,mu_comp
INTEGER :: i,j,k,LJ_N
IF (.NOT. PRESENT(lbounds)) THEN
lb = -3_MK
ELSE
lb = lbounds
END IF
IF (.NOT. PRESENT(ubounds)) THEN
ub = 3_MK
ELSE
ub = ubounds
END IF
!-------------------------------------------------------------------------
! Lennard Jones Parameters
!-------------------------------------------------------------------------
! Lennard-Jones Potential
!-------------------------------------------------------------------------
LJsigma=1.0_MK
LJeps=1.0_MK
! Compression weight
mu_comp = options%LJ_comp*LJeps
! number of particles
LJ_N = (m+6)/3
res(:) = 0.0_MK
! loop over all particle configurations
! Sequential computation of n solutions
DO i=1,n
! place first particle at the origin
part_mat(:,:)=0.0_MK
! place second particle at the x-axis
part_mat(2,1)=vars(1,i)
part_mat(2,2)=0.0_MK
part_mat(2,3)=0.0_MK
! place third particle in the x-y plane
part_mat(3,1)=vars(2,i)
part_mat(3,2)=vars(3,i)
part_mat(3,3)=0.0_MK
! build particle position matrix
DO j=4,m-2,3
k=(j+8)/3
part_mat(k,1)=vars(j,i)
part_mat(k,2)=vars(j+1,i)
part_mat(k,3)=vars(j+2,i)
END DO
! initialize distance
dist=0.0_MK
! compute distance matrix and the Lennard-Jones energy
DO k=1,LJ_N-1
DO l=k+1,LJ_N
dist=(sum((part_mat(k,:)-part_mat(l,:))**2))**(0.5)
res(i) = res(i) + 4.0_MK*LJeps*((LJsigma/dist)**12 - (LJsigma/dist)**6)
END DO
END DO
part_mean(1,1) = 1./LJ_N*sum(part_mat(:,1))
part_mean(1,2) = 1./LJ_N*sum(part_mat(:,2))
part_mean(1,3) = 1./LJ_N*sum(part_mat(:,3))
DO k=1,LJ_N
res(i) = res(i) + mu_comp * sum((part_mat(k,:)-part_mean(1,:))**2) / LJsigma**2
END DO
END DO
END SUBROUTINE LJ_POT_COMP
SUBROUTINE LJ_POT(res,vars,m,n,lbounds,ubounds)
USE cmaes_param_mod
!USE deriv_class ! make the module accessible
!Parameters
REAL(MK), DIMENSION(n), INTENT(out) :: res
REAL(MK), DIMENSION(m,n),INTENT(in) :: vars
INTEGER, INTENT(in) :: m ! (3*LJ_N-6) dim vector
INTEGER, INTENT(in) :: n
REAL(MK), OPTIONAL :: lbounds
REAL(MK), OPTIONAL :: ubounds
! Locally used variables
REAL(MK) :: dist
REAL(MK), DIMENSION((m+6)/3,3) :: part_mat
REAL(MK), DIMENSION(1,3) :: part_mean
REAL(MK) :: ub,lb
REAL(MK) :: LJsigma,LJeps
INTEGER :: i,j,k,LJ_N
IF (.NOT. PRESENT(lbounds)) THEN
lb = -3_MK
ELSE
lb = lbounds
END IF
IF (.NOT. PRESENT(ubounds)) THEN
ub = 3_MK
ELSE
ub = ubounds
END IF
!-------------------------------------------------------------------------
! Lennard Jones Parameters
!-------------------------------------------------------------------------
! Lennard-Jones Potential
!-------------------------------------------------------------------------
LJsigma=1.0_MK
LJeps=1.0_MK
! number of particles
LJ_N = (m+6)/3
res(:) = 0.0_MK
! loop over all particle configurations
! Sequential computation of n solutions
DO i=1,n
! place first particle at the origin
part_mat(:,:)=0.0_MK
! place second particle at the x-axis
part_mat(2,1)=vars(1,i)
part_mat(2,2)=0.0_MK
part_mat(2,3)=0.0_MK
! place third particle in the x-y plane
part_mat(3,1)=vars(2,i)
part_mat(3,2)=vars(3,i)
part_mat(3,3)=0.0_MK
! build particle position matrix
DO j=4,m-2,3
k=(j+8)/3
part_mat(k,1)=vars(j,i)
part_mat(k,2)=vars(j+1,i)
part_mat(k,3)=vars(j+2,i)
END DO
! initialize distance
dist=0.0_MK
! compute distance matrix and the Lennard-Jones energy
DO k=1,LJ_N-1
DO l=k+1,LJ_N
dist=(sum((part_mat(k,:)-part_mat(l,:))**2))**(0.5)
res(i) = res(i) + 4.0_MK*LJeps*((LJsigma/dist)**12 - (LJsigma/dist)**6)
END DO
END DO
END DO
END SUBROUTINE LJ_POT
This diff is collapsed.
SUBROUTINE LJ_write_out(filename,X,m,n) !,V,ELJ,GTEST,SECT)
USE cmaes_param_mod
USE cmaes_mod
IMPLICIT NONE
CHARACTER(LEN = 200),INTENT(IN) :: filename
REAL(MK), DIMENSION(m,n),INTENT(in) :: X
INTEGER, INTENT(in) :: m ! (3*LJ_N-6) dim vector
INTEGER, INTENT(in) :: n
!local vars
REAL(MK), DIMENSION(m+6) :: X_full
REAL(MK), DIMENSION(m+6) :: V_full
INTEGER :: i
LOGICAL :: GTEST
! local variables
INTEGER,DIMENSION(:),ALLOCATABLE :: Npart_vec
REAL(MK),DIMENSION(:,:),ALLOCATABLE :: xp_iproc
INTEGER :: ip,count,iproc
CHARACTER :: altloc
CHARACTER (LEN = 4 ) :: atom_name
CHARACTER :: chains
CHARACTER (LEN = 2 ) :: charge
CHARACTER (LEN = 2 ) :: element
CHARACTER :: icode
INTEGER :: resno
REAL(8) :: occ
CHARACTER (LEN = 3 ) :: resname
CHARACTER (LEN = 4 ) :: segid
INTEGER :: info
!====================================================================!
! initialize
atom_name = 'H '
altloc = ' '
resname = ' '
chains = ' '
resno = 0
icode = ' '
occ = 1.0D+00
segid = ' '
element = ' '
charge = ' '
X_full = 0
X_full(4) = X(1,1)
X_full(7) = X(2,1)
X_full(8) = X(3,1)
DO i = 4,m
X_full(i+6) = X(i,1)
END DO
X_full(4) = tmptest(1)
X_full(7) = tmptest(2)
X_full(8) =tmptest(3)
DO i = 4,m
X_full(i+6) = tmptest(i)
!X_full(i+6) = X(i,n)
END DO
!ATOM 1607 N LYS A 233 -26.180 4.759 -18.385 1.00 41.53 N
!====================================================================!
! open file
IF (countIter .EQ. 1) THEN
OPEN(99,FILE=filename,FORM='FORMATTED',STATUS='REPLACE',IOSTAT=info)
IF (info .NE. 0) THEN
WRITE(*,*)'ERROR in lj_writeout: opening file failed.'
GOTO 9999
ENDIF
END IF
WRITE (99,'(a6,i5)')'MODEL ',countRestarts
DO i = 1,m+6,3
WRITE (99, &
'(a6,i5,1x,a4,a1,a3,a1,1x,i4,a1,3x,3f8.3,2f6.2,6x,a4,a2,a2)' ) &
'ATOM ', i, atom_name, altloc, resname, chains, &
& resno,icode,X_full(i),X_full(i+1),X_full(i+2),occ,0d0,segid, &
& element, charge
END DO
!====================================================================!
! terminate pdb-file and close file
WRITE (99,'(a6)')'ENDMDL'
!WRITE (99, &
! '(a6,i5,1x,4x,a1,a3,1x,a1,i4,a1)' ) &
! 'TER ', ip+1, altloc, resname, chains, resno, icode
!CLOSE(99,IOSTAT=info)
!IF (info .NE. 0) THEN
! WRITE(*,*)'ERROR in lj_writeout: closing file failed.'
! GOTO 9999
!ENDIF
9999 CONTINUE ! jump here upon error
END SUBROUTINE LJ_write_out
This diff is collapsed.
SUBROUTINE water_write_out(filename,X,m,n) !,V,ELJ,GTEST,SECT)
USE cmaes_param_mod
USE cmaes_mod
IMPLICIT NONE
CHARACTER(LEN = 200),INTENT(IN) :: filename
REAL(MK), DIMENSION(m,n),INTENT(in) :: X
INTEGER, INTENT(in) :: m ! (3*LJ_N-6) dim vector
INTEGER, INTENT(in) :: n
!local vars
REAL(MK), DIMENSION(m/2*3) :: X_full
!REAL(MK), DIMENSION(m+6) :: V_full
INTEGER :: i,j
LOGICAL :: GTEST
! local variables
INTEGER,DIMENSION(:),ALLOCATABLE :: Npart_vec
REAL(MK),DIMENSION(:,:),ALLOCATABLE :: xp_iproc
INTEGER :: ip,count,iproc
CHARACTER :: altloc
CHARACTER (LEN = 4 ) :: atom_name
CHARACTER :: chains
CHARACTER (LEN = 2 ) :: charge
CHARACTER (LEN = 2 ) :: element
CHARACTER :: icode
INTEGER :: resno
REAL(8) :: occ
CHARACTER (LEN = 3 ) :: resname
CHARACTER (LEN = 4 ) :: segid
INTEGER :: info
REAL(MK),DIMENSION(9) :: Cords
!====================================================================!
! initialize
altloc = ' '
resname = ' '
chains = ' '
resno = 0
icode = ' '
occ = 1.0D+00
segid = ' '
element = ' '
charge = ' '
DO i = 1,m/2,3
CALL TIPIO( X(i,1), X(i+1,1), X(i+2,1),&
X(m/2+i,1),X(m/2+1+i,1),X(m/2+2+i,1), Cords)
DO j = 1,9
X_full((i-1)/3*9+j) = Cords(j)
END DO
END DO
!ATOM 1607 N LYS A 233 -26.180 4.759 -18.385 1.00 41.53 N
!====================================================================!
! open file
IF (countIter .EQ. 1) THEN
OPEN(99,FILE=filename,FORM='FORMATTED',STATUS='REPLACE',IOSTAT=info)
IF (info .NE. 0) THEN
WRITE(*,*)'ERROR in water_writeout: opening file failed.'
GOTO 9999
ENDIF
END IF
WRITE (99,'(a6,i5)')'MODEL ',countIter + 1
DO i = 1,m/2*3,9
atom_name = 'O '
WRITE (99, &
'(a6,i5,1x,a4,a1,a3,a1,1x,i4,a1,3x,3f8.3,2f6.2,6x,a4,a2,a2)' ) &
'ATOM ', i, atom_name, altloc, resname, chains, &
& resno,icode,X_full(i),X_full(i+1),X_full(i+2),occ,0d0,segid, &
& element, charge
atom_name = 'H '
WRITE (99, &
'(a6,i5,1x,a4,a1,a3,a1,1x,i4,a1,3x,3f8.3,2f6.2,6x,a4,a2,a2)' ) &
'ATOM ', i+1, atom_name, altloc, resname, chains, &
& resno,icode,X_full(i+3),X_full(i+4),X_full(i+5),occ,0d0,segid, &
& element, charge
WRITE (99, &
'(a6,i5,1x,a4,a1,a3,a1,1x,i4,a1,3x,3f8.3,2f6.2,6x,a4,a2,a2)' ) &
'ATOM ', i+2, atom_name, altloc, resname, chains, &
& resno,icode,X_full(i+6),X_full(i+7),X_full(i+8),occ,0d0,segid, &
& element, charge
END DO
!====================================================================!
! terminate pdb-file and close file
WRITE (99,'(a6)')'ENDMDL'
!WRITE (99, &
! '(a6,i5,1x,4x,a1,a3,1x,a1,i4,a1)' ) &
! 'TER ', ip+1, altloc, resname, chains, resno, icode
!CLOSE(99,IOSTAT=info)
!IF (info .NE. 0) THEN
! WRITE(*,*)'ERROR in lj_writeout: closing file failed.'
! GOTO 9999
!ENDIF
9999 CONTINUE ! jump here upon error
END SUBROUTINE water_write_out
SUBROUTINE TIPIO(X1, Y1, Z1, L1, M1, N1, COORDS)
IMPLICIT NONE
DOUBLE PRECISION X1, Y1, Z1, COORDS(*), C2A1, M1, L1, N1, ALPHA1, CA1, S1, C3A1, L12, M12, N12
L12=L1**2
M12=M1**2
N12=N1**2
ALPHA1=SQRT(L12+M12+N12)
CA1=COS(ALPHA1)
C2A1=CA1
IF (ALPHA1.LT.0.0001D0) THEN
! C3A1=(-ALPHA1/2+ALPHA1**3/24)
C3A1=(-0.5D0+ALPHA1**2/24.0D0)
S1=(1.0D0-ALPHA1**2/6)
ELSE
C3A1=(CA1-1.0D0)/ALPHA1**2
S1=SIN(ALPHA1)/ALPHA1
ENDIF
COORDS(1) = X1
COORDS(2) = Y1
COORDS(3) = Z1
COORDS(4) = 0.756950327*c2a1 - c3a1*l1*(0.756950327*l1 - 0.585882276*n1) + 0.585882276*m1*s1 + X1
COORDS(5) = -(c3a1*m1*(0.756950327*l1 - 0.585882276*n1)) + (-0.585882276*l1 - 0.756950327*n1)*s1 + Y1
COORDS(6) = -0.585882276*c2a1 - c3a1*(0.756950327*l1 - 0.585882276*n1)*n1 + 0.756950327*m1*s1 + Z1
COORDS(7) = -0.756950327*c2a1 + c3a1*l1*(0.756950327*l1 + 0.585882276*n1) + 0.585882276*m1*s1 + X1
COORDS(8) = c3a1*m1*(0.756950327*l1 + 0.585882276*n1) + (-0.585882276*l1 + 0.756950327*n1)*s1 + Y1
COORDS(9) = -0.585882276*c2a1 + c3a1*(0.756950327*l1 + 0.585882276*n1)*n1 - 0.756950327*m1*s1 + Z1
RETURN
END
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment