[CPMD-list] constraints revisited
Axel Kohlmeyer
axel.kohlmeyer at theochem.ruhr-uni-bochum.de
Sat Nov 2 18:54:19 CET 2002
On Sat, 02 Nov 2002 18:18:05 +0100 Phineus Markwick wrote:
> This is a multi-part message in MIME format.
> --------------010009080001030009090302
> Content-Type: text/plain; charset=us-ascii; format=flowed
> Content-Transfer-Encoding: 7bit
>
> Dear cpmd-list,
hi phineus!
>
> In trying to incorporate a new constraint into cpmd, I have a problem in
> the program CNSTFC.F at the point:
>
> DO K=1,KMAX
> DO N=1,NODIM
> ANORM(N,I)=ANORM(N,I)+DX(K)*ASKEL(N,I,K)
> ENDDO
> ENDDO
> ENDDO
>
> can anybody tell me exactly what the matrix ASKEL(N,I,K) is and where
> its components are defined?
hmm, i really don't know a lot about cpmd internals, but this question
can hopefully be answered anyway:
[25|18:45]~/compile/cpmd/SOURCE> grep -n -i askel *.F *.inc
cnstfc.F:92: ANORM(N,I)=ANORM(N,I)+DX(K)*ASKEL(N,I,K)
detdof.F:114: CALL MEMORY(IP_ASKEL,12*NODIM*MCNSTR,'ASKEL')
detdof.F:121: CALL AZZERO(ASKEL,12*NODIM*MCNSTR)
detdof.F:228: IF(JX.NE.0) ASKEL(JX,I,NN+1)=1.0D0
detdof.F:229: IF(JY.NE.0) ASKEL(JY,I,NN+2)=1.0D0
detdof.F:230: IF(JZ.NE.0) ASKEL(JZ,I,NN+3)=1.0D0
detdof.F:245: IF(JX.NE.0) ASKEL(JX,I,NN+1)=FF
detdof.F:246: IF(JY.NE.0) ASKEL(JY,I,NN+2)=FF
detdof.F:247: IF(JZ.NE.0) ASKEL(JZ,I,NN+3)=FF
cotr.inc:44:C == ASKEL(NODIM,MCNSTR,12)
==
cotr.inc:56: & CNSVAL(*),ANORM(NODIM,*),ASKEL(NODIM,MCNSTR,*),
cotr.inc:62: * (IP_ANORM,ANORM),(IP_ASKEL,ASKEL),(IP_FCSTR,FCSTR
),
cotr.inc:68: * IP_NTCNST,IP_CNSVAL,IP_ANORM,IP_ASKEL,IP_FCSTR,
>
> As I am now sure that the problem is in the program CNSTFC.F (thanks
> Axel!), I have included my edited version of this program as an attachment.
> In my test system, the distance constraint is formulated using the
> coordinates of 10 atoms, so dX (or dR) is a vector with 30 components,
> and KMAX=30. I think that the problem comes from the fact that all
> N*I*30 components of the matrix ASKEL are not defined.
>
> Once again, I would appreciate any help or comments,
now let's see what you have changed:
[26|18:46]~/compile/cpmd/SOURCE> diff -u cnstfc.F ~/work/cnstfc.F
--- cnstfc.F Tue Nov 14 13:21:19 2000
+++ /home/akohlmey/work/cnstfc.F Sat Nov 2 18:42:22 2002
@@ -14,7 +14,9 @@
LOGICAL LAGRA
C Variables
INTEGER ISUB,IA,IB,IC,ID,ITYP,I,J,K,N,KMAX,IX,IDAMAX
- REAL*8 X1(3),X2(3),X3(3),X4(3),DX(12),CVAL
+ INTEGER IE,IFF,IG,IH,II,IJ
+ REAL*8 X1(3),X2(3),X3(3),X4(3),DX(30),CVAL
+ REAL*8 X5(3),X6(3),X7(3),X8(3),X9(3),X10(3)
hmm, the DX(12) seems to correspond with the definition of askel in
detdof.F:114. but let's continue:
@@ -60,6 +68,23 @@
CALL FUNCD(FC(I),FV(I),CVAL,X1,X2)
CALL DIFFD(DX,X1,X2)
KMAX=6
+ ELSEIF(ITYP.EQ.9) THEN
+C..NEW CONSTRAINT
+ CALL FILLC(IA,TSCR,X1)
+ CALL FILLC(IB,TSCR,X2)
+ CALL FILLC(IC,TSCR,X3)
+ CALL FILLC(ID,TSCR,X4)
+ CALL FILLC(IE,TSCR,X5)
+ CALL FILLC(IFF,TSCR,X6)
+ CALL FILLC(IG,TSCR,X7)
+ CALL FILLC(IH,TSCR,X8)
+ CALL FILLC(II,TSCR,X9)
+ CALL FILLC(IJ,TSCR,X10)
+ CALL FUNCCONS(FC(I),FV(I),CVAL,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10)
+ CALL DIFFCONS(DX,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10)
+C PRINT*,'DX=',DX
+C PRINT*,'XS=',X1,X2,X3,X4,X5,X6,X7,X8,X9,X10
+ KMAX=30
ELSEIF(ITYP.EQ.5) THEN
ok, here you change kmax from 6 to 30, but the comment in cotr.inc:44
seems to indicate a limit of 12. oops. i would say, try replacing the
in 12 detdof.F:114 and detdof.F:121 with 30. but someone who _really_
knows cpmd should confirm that.
cheers,
axel.
>
> best regards,
> Phineus Markwick.
>
> --------------010009080001030009090302
> Content-Type: text/plain;
> name="cnstfc.F"
> Content-Transfer-Encoding: 7bit
> Content-Disposition: inline;
> filename="cnstfc.F"
>
> C ==================================================================
> SUBROUTINE CNSTFC(TAU0,TSCR,DXPAR,CNMAX,LAGRA)
> C ==--------------------------------------------------------------==
> IMPLICIT NONE
> INCLUDE 'system.h'
> INCLUDE 'ions.inc'
> INCLUDE 'cotr.inc'
> INCLUDE 'rmas.inc'
> INCLUDE 'tpar.inc'
> INCLUDE 'adat.inc'
> INCLUDE 'ener.inc'
> C Arguments
> REAL*8 TAU0(3,NAX,NSX),TSCR(3,NAX*NSX),DXPAR(*),CNMAX
> LOGICAL LAGRA
> C Variables
> INTEGER ISUB,IA,IB,IC,ID,ITYP,I,J,K,N,KMAX,IX,IDAMAX
> INTEGER IE,IFF,IG,IH,II,IJ
> REAL*8 X1(3),X2(3),X3(3),X4(3),DX(30),CVAL
> REAL*8 X5(3),X6(3),X7(3),X8(3),X9(3),X10(3)
> C ==--------------------------------------------------------------==
> CALL TISET(' CNSTFC',ISUB)
> CNMAX=0.0D0
> C No constraints -> return.
> IF(MCNSTR.EQ.0) RETURN
> CALL DUM2(TAU0,TSCR)
> CALL AZZERO(ANORM,MCNSTR*NODIM)
> DO I=1,MCNSTR
> CVAL=CNSVAL(I)
> ITYP=NTCNST(1,I)
> IA=NTCNST(2,I)
> IB=NTCNST(3,I)
> IC=NTCNST(4,I)
> ID=NTCNST(5,I)
> IE=NTCNST(6,I)
> IFF=NTCNST(7,I)
> IG=NTCNST(8,I)
> IH=NTCNST(9,I)
> II=NTCNST(10,I)
> IJ=NTCNST(11,I)
> IF(ITYP.EQ.1) THEN
> C..stretch
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FUNCR(FC(I),FV(I),CVAL,X1,X2)
> CALL DIFFR(DX,X1,X2)
> KMAX=6
> C..angle
> ELSEIF(ITYP.EQ.2) THEN
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FILLC(IC,TSCR,X3)
> CALL FUNCT(FC(I),FV(I),CVAL,X1,X2,X3)
> CALL DIFFT(DX,X1,X2,X3)
> KMAX=9
> ELSEIF(ITYP.EQ.3) THEN
> C..dihedral angle
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FILLC(IC,TSCR,X3)
> CALL FILLC(ID,TSCR,X4)
> CALL FUNCO(FC(I),FV(I),CVAL,X1,X2,X3,X4)
> CALL DIFFO(DX,X1,X2,X3,X4)
> KMAX=12
> ELSEIF(ITYP.EQ.4) THEN
> C..distance
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FUNCD(FC(I),FV(I),CVAL,X1,X2)
> CALL DIFFD(DX,X1,X2)
> KMAX=6
> ELSEIF(ITYP.EQ.9) THEN
> C..NEW CONSTRAINT
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FILLC(IC,TSCR,X3)
> CALL FILLC(ID,TSCR,X4)
> CALL FILLC(IE,TSCR,X5)
> CALL FILLC(IFF,TSCR,X6)
> CALL FILLC(IG,TSCR,X7)
> CALL FILLC(IH,TSCR,X8)
> CALL FILLC(II,TSCR,X9)
> CALL FILLC(IJ,TSCR,X10)
> CALL FUNCCONS(FC(I),FV(I),CVAL,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10)
> CALL DIFFCONS(DX,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10)
> C PRINT*,'DX=',DX
> C PRINT*,'XS=',X1,X2,X3,X4,X5,X6,X7,X8,X9,X10
> KMAX=30
> ELSEIF(ITYP.EQ.5) THEN
> C..out of plane
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FILLC(IC,TSCR,X3)
> CALL FILLC(ID,TSCR,X4)
> CALL FUNCP(FC(I),FV(I),CVAL,X1,X2,X3,X4)
> CALL DIFFP(DX,X1,X2,X3,X4)
> KMAX=12
> ELSEIF(ITYP.EQ.6) THEN
> C..coordination number
> CALL COORNUM(IA,TSCR,NAT,ANORM(1,I),LSKPTR,
> * FC(I),FV(I),CVAL)
> KMAX=0
> ELSEIF(ITYP.EQ.7) THEN
> C..difference between distances- 3 atoms
> CALL FILLC(IA,TSCR,X1)
> CALL FILLC(IB,TSCR,X2)
> CALL FILLC(IC,TSCR,X3)
> CALL FUNCDD(FC(I),FV(I),CVAL,X1,X2,X3)
> CALL DIFFDD(DX,X1,X2,X3)
> KMAX=9
> ELSE
> WRITE(*,*) ' CNSTFC! I=',I,' TYPE=',ITYP
> CALL STOPGM('CNSTFC','UNKNOWN TYPE OF CONSTRAINT')
> END IF
>
> PRINT*,'hallo world 1'
>
> DO K=1,KMAX
> DO N=1,NODIM
> ANORM(N,I)=ANORM(N,I)+DX(K)*ASKEL(N,I,K)
> ENDDO
> ENDDO
> ENDDO
>
> PRINT*,'hallo world 2'
>
> ECNSTR=0.0D0
> CALL AZZERO(FCSTR,NODIM)
> DO I=1,MCNSTR
> DO J=1,NODIM
> IF(LAGRA) THEN
> FCSTR(J)=FCSTR(J)+XLAGR(I)*ANORM(J,I)
> ELSE
> FCSTR(J)=FCSTR(J)+CSIGM(I)*ANORM(J,I)*FC(I)
> ECNSTR=ECNSTR+0.5D0*CSIGM(I)*FC(I)*FC(I)
> ENDIF
> ENDDO
> ENDDO
> IX=IDAMAX(NODIM,FCSTR(1),1)
> CNMAX=ABS(FCSTR(IX))
> CALL DAXPY(NODIM,1.0D0,FCSTR(1),1,DXPAR(1),1)
> CALL TIHALT(' CNSTFC',ISUB)
> C ==--------------------------------------------------------------==
> RETURN
> END
> C ==================================================================
> SUBROUTINE FILLC(IAT,TAU,X)
> C ==--------------------------------------------------------------==
> C == Extract the coordinates X(1:3) of IAT index in TAU (TSCR) ==
> C == Works also if dummy atoms ==
> C ==--------------------------------------------------------------==
> IMPLICIT NONE
> INCLUDE 'system.h'
> INCLUDE 'ions.inc'
> INCLUDE 'cotr.inc'
> C Arguments
> INTEGER IAT
> REAL*8 TAU(3,NAX*NSX),X(3)
> C Local
> INTEGER ID,NAA,ITYP
> C ==--------------------------------------------------------------==
> IF(IAT.LE.0) THEN
> C Do nothing
> RETURN
> ELSEIF(IAT.LE.NAT) THEN
> C Real atoms
> X(1)=TAU(1,IAT)
> X(2)=TAU(2,IAT)
> X(3)=TAU(3,IAT)
> ELSEIF(IAT.LE.NAT+NDAT) THEN
> C Dummy atoms (type 2).
> NAA=IAT-NAT
> ITYP=LISTDA(NAA,1)
> ID=LISTDA(NAA,2)
> IF(ITYP.EQ.2) THEN
> X(1)=DUMMY2(1,ID)
> X(2)=DUMMY2(2,ID)
> X(3)=DUMMY2(3,ID)
> ENDIF
> ENDIF
> C ==--------------------------------------------------------------==
> RETURN
> END
> C ==================================================================
> SUBROUTINE COORNUM(IA,TSCR,NAT,AN,LSK,FCI,FVI,CVAL)
> C ==--------------------------------------------------------------==
> IMPLICIT NONE
> INCLUDE 'system.h'
> INCLUDE 'cotr.inc'
> C Arguments
> INTEGER IA,NAT,LSK(3,NAT)
> REAL*8 TSCR(3,*),AN(*),FCI,FVI,CVAL
> C Variables
> REAL*8 X0,Y0,Z0,DX,DY,DZ,DD,FF,DF
> INTEGER IAT,K,L1,L2,L3
> C ==--------------------------------------------------------------==
> L1 = LSK(1,IA)
> L2 = LSK(2,IA)
> L3 = LSK(3,IA)
> X0 = TSCR(1,IA)
> Y0 = TSCR(2,IA)
> Z0 = TSCR(3,IA)
> FVI = 0.D0
> DO IAT=1,NAT
> DX=TSCR(1,IAT)-X0
> DY=TSCR(2,IAT)-Y0
> DZ=TSCR(3,IAT)-Z0
> CALL PBC(DX,DY,DZ,DX,DY,DZ,1,APBC,IBRAV)
> DD=SQRT(DX*DX+DY*DY+DZ*DZ)
> C..exclude self interaction
> IF(DD.GT.1.D-2) THEN
> DF=C_KAPPA*(DD-C_RC)
> FVI=FVI+1.D0/(EXP(DF)+1.D0)
> FF=-0.5*C_KAPPA/(COSH(DF)+1.D0)/DD
> IF(LSK(1,IAT).NE.0) THEN
> K=LSK(1,IAT)
> AN(K)=AN(K)+FF*DX
> ENDIF
> IF(LSK(2,IAT).NE.0) THEN
> K=LSK(2,IAT)
> AN(K)=AN(K)+FF*DY
> ENDIF
> IF(LSK(3,IAT).NE.0) THEN
> K=LSK(3,IAT)
> AN(K)=AN(K)+FF*DZ
> ENDIF
> IF(L1.NE.0) AN(L1)=AN(L1)-FF*DX
> IF(L2.NE.0) AN(L2)=AN(L2)-FF*DY
> IF(L3.NE.0) AN(L3)=AN(L3)-FF*DZ
> ENDIF
> ENDDO
> FCI=FVI-CVAL
> C ==--------------------------------------------------------------==
> RETURN
> END
> C ==================================================================
>
> --------------010009080001030009090302--
>
> _______________________________________________
> CPMD-list mailing list
> CPMD-list at cpmd.org
> http://www.cpmd.org/mailman/listinfo/cpmd-list
>
--
=======================================================================
Axel Kohlmeyer e-mail: axel.kohlmeyer at theochem.ruhr-uni-bochum.de
Lehrstuhl fuer Theoretische Chemie Phone: ++49 (0)234/32-26673
Ruhr-Universitaet Bochum - NC 03/53 Fax: ++49 (0)234/32-14045
D-44780 Bochum http://www.theochem.ruhr-uni-bochum.de
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.
More information about the CPMD-list
mailing list