5 1.00 1.500 2.000 ! Dbz, HD
2 ! Aib, HB
3 2.090 ! Abu, HG
- 5 1.00 1.500 2.000 ! Dbz, HD
- 2 ! Aib, HB
- 3 2.090 ! Abu, HG
--- /dev/null
+19
+************ p-cys
+ 4
+ 0 20.9855397362
+ 1 23.2362754541
+ 2 15.2145213599
+ 3 5.8457062573
+ 4 2.2247834920
+************ cys-p
+ 4
+ 0 30.5588852886
+ 1 38.9616272470
+ 2 26.2903926480
+ 3 10.9408066120
+ 4 4.1305862378
+************ p-met
+ 4
+ 0 14.6528703204
+ 1 10.8322908087
+ 2 6.4948703253
+ 3 2.5156972424
+ 4 0.7678847831
+************ met-p
+ 4
+ 0 23.7090964872
+ 1 26.8032006361
+ 2 18.7519376751
+ 3 9.2894391762
+ 4 3.2385892828
+************ p-phe
+ 4
+ 0 13.0676515444
+ 1 9.1621853821
+ 2 5.8535312419
+ 3 2.6234191522
+ 4 0.5272694493
+************ phe-p
+ 4
+ 0 19.7887293956
+ 1 20.6087917953
+ 2 14.5650743251
+ 3 8.2456178099
+ 4 2.6755969113
+************ p-ile
+ 4
+ 0 21.6748065898
+ 1 22.4712439569
+ 2 12.3323689766
+ 3 3.7675300327
+ 4 0.8867786667
+************ ile-p
+ 4
+ 0 53.8375910531
+ 1 77.1108545106
+ 2 51.4737848309
+ 3 22.2224900972
+ 4 7.0592435547
+************ p-leu
+ 4
+ 0 26.2693630015
+ 1 31.4777831332
+ 2 20.2250151158
+ 3 8.8864378397
+ 4 2.9300820063
+************ leu-p
+ 4
+ 0 39.6311445722
+ 1 54.4209827652
+ 2 35.3763863718
+ 3 15.5966073687
+ 4 4.4475409968
+************ p-val
+ 4
+ 0 17.6284319040
+ 1 17.5807828769
+ 2 8.6572369364
+ 3 2.0219350588
+ 4 0.3680322524
+************ val-p
+ 4
+ 0 67.1316515604
+ 1 101.2634978135
+ 2 69.9320403365
+ 3 30.3867583774
+ 4 10.6357698795
+************ p-trp
+ 4
+ 0 11.6095328599
+ 1 6.0406520753
+ 2 4.2836365803
+ 3 2.2395796734
+ 4 0.6397809424
+************ trp-p
+ 4
+ 0 19.7869247023
+ 1 19.0462135491
+ 2 14.5261229711
+ 3 8.9346155113
+ 4 3.3332821409
+************ p-tyr
+ 4
+ 0 10.6110403221
+ 1 5.1063715761
+ 2 2.9094288206
+ 3 1.2858200387
+ 4 0.0650991057
+************ tyr-p
+ 4
+ 0 16.3834307203
+ 1 14.5187824862
+ 2 10.2204383636
+ 3 6.0436432257
+ 4 1.8640778923
+************ p-ala
+ 4
+ 0 25.5529476036
+ 1 26.8929521099
+ 2 26.4141845455
+ 3 8.6066730099
+ 4 8.1703930067
+************ ala-p
+ 4
+ 0 158.5960056306
+ 1 260.5554552766
+ 2 173.6920621211
+ 3 73.5801296687
+ 4 21.9400038548
+************ p-thr
+ 4
+ 0 14.2192826135
+ 1 8.7197768861
+ 2 3.7577261897
+ 3 -1.0782876908
+ 4 -1.1274709648
+************ thr-p
+ 4
+ 0 45.0221273479
+ 1 60.6096030190
+ 2 43.4909169091
+ 3 17.9626982908
+ 4 6.6133607156
+************ p-ser
+ 4
+ 0 13.5839488564
+ 1 6.6706098352
+ 2 4.6534430546
+ 3 0.2492254367
+ 4 0.0018838653
+************ ser-p
+ 4
+ 0 53.3414876859
+ 1 78.7194975176
+ 2 56.0713187851
+ 3 26.1372171643
+ 4 9.2684992822
+************ p-gln
+ 4
+ 0 17.0265171919
+ 1 12.2729170319
+ 2 8.4659874420
+ 3 3.3060953614
+ 4 0.7613430027
+************ gln-p
+ 4
+ 0 18.9455470339
+ 1 16.1932316164
+ 2 11.8150807911
+ 3 5.5714895281
+ 4 1.9069932595
+************ p-asn
+ 4
+ 0 15.0667382548
+ 1 9.0217864886
+ 2 5.9497922773
+ 3 1.2583204641
+ 4 0.6835557859
+************ asn-p
+ 4
+ 0 30.4827979254
+ 1 37.0561579432
+ 2 25.5895542308
+ 3 11.9135799219
+ 4 3.7539802438
+************ p-glu
+ 4
+ 0 15.6056591126
+ 1 11.4184667608
+ 2 7.6406356194
+ 3 3.1869364269
+ 4 0.8055606024
+************ glu-p
+ 4
+ 0 17.4767844966
+ 1 14.1992440205
+ 2 10.2588118394
+ 3 5.1122037666
+ 4 1.6769129114
+************ p-asp
+ 4
+ 0 13.8330918514
+ 1 8.2889193990
+ 2 5.3674424779
+ 3 1.4951210044
+ 4 0.6724647127
+************ asp-p
+ 4
+ 0 31.3411709322
+ 1 38.5241372893
+ 2 26.4391050476
+ 3 12.7795122491
+ 4 3.8993221861
+************ p-his
+ 4
+ 0 15.3664120018
+ 1 10.9142271706
+ 2 7.0027777295
+ 3 2.6145603184
+ 4 0.3733345459
+************ his-p
+ 4
+ 0 19.9958016121
+ 1 19.9841853554
+ 2 14.2411214420
+ 3 7.3269627286
+ 4 2.2235893664
+************ p-arg
+ 4
+ 0 15.6841789437
+ 1 8.5094056867
+ 2 7.1651416581
+ 3 3.7459223928
+ 4 0.7401732041
+************ arg-p
+ 4
+ 0 17.6421139033
+ 1 10.9172352462
+ 2 9.4536664737
+ 3 5.4955136323
+ 4 1.9054031555
+************ p-lys
+ 4
+ 0 16.8219722502
+ 1 14.1766716664
+ 2 9.2479297479
+ 3 4.1072620109
+ 4 1.1847542958
+************ lys-p
+ 4
+ 0 17.6927902003
+ 1 15.9785187922
+ 2 11.0246553123
+ 3 5.4292106221
+ 4 1.5303688445
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+************ p-pro
+ 4
+ 0 241.2416350955
+ 1 -96.3545930066
+ 2 284.1188810874
+ 3 -12.8583927986
+ 4 49.9760710382
+************ pro-p
+ 4
+ 0 360.3641334040
+ 1 598.3562129486
+ 2 409.3227118177
+ 3 175.4335390823
+ 4 55.9961827971
+
endif(UNRES_MD_FF STREQUAL "GAB")
if(UNRES_NEWGRAD)
- set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG ")
+# set(CPPFLAGS "${CPPFLAGS} -DFIVEDIAG ")
+set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG -DSC_END")
endif()
+# set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG ")
#=========================================
!c-----------------------------------------------------------------
!c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups
!c inside, implemented with five-diagonal inertia matrix
+ use energy_data
implicit none
- include 'DIMENSIONS'
- include 'COMMON.VAR'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.MD'
- include 'COMMON.LAGRANGE.5diag'
- include 'COMMON.IOUNITS'
- double precision KE_total
- integer i,j,k,iti
- double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
- mag1,mag2,v(3)
+ real(kind=8):: KE_total,KEt_p,KEt_sc,KEr_p,KEr_sc,mag1,mag2
+ integer i,j,k,iti,mnum
+ real(kind=8),dimension(3) :: incr,v
KEt_p=0.0d0
KEt_sc=0.0d0
incr(j)=d_t(j,0)
enddo
do i=nnt,nct-1
+ mnum=molnum(i)
!c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
!c Skip dummy peptide groups
- if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
+ if (itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
do j=1,3
v(j)=incr(j)+0.5d0*d_t(j,i)
enddo
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
!c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
- KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
endif
do j=1,3
incr(j)=incr(j)+d_t(j,i)
incr(j)=d_t(j,0)
enddo
do i=nnt,nct
- iti=iabs(itype(i))
- if (itype(i).eq.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or.mnum.ge.3) then
do j=1,3
v(j)=incr(j)
enddo
v(j)=incr(j)+d_t(j,nres+i)
enddo
endif
-!c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
-!c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
- KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+ write (iout,*) "i",i," msc",msc(iti,1)," v",(v(j),j=1,3)
+ KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
do j=1,3
incr(j)=incr(j)+d_t(j,i)
! write(iout,*) 'KEt_sc', KEt_sc
! The part due to stretching and rotation of the peptide groups
do i=nnt,nct-1
- if (itype(i).ne.ntyp1.and.itype(i+1).ne.ntyp1) then
+ mnum=molnum(i)
+ if (itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
+ if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
! write (iout,*) "i",i
! write (iout,*) "i",i," mag1",mag1," mag2",mag2
do j=1,3
incr(j)=d_t(j,i)
enddo
!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
- KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
+ KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
+incr(3)*incr(3))
endif
enddo
!c write(iout,*) 'KEr_p', KEr_p
!c The rotational part of the side chain virtual bond
do i=nnt,nct
- iti=iabs(itype(i))
- if (itype(i).ne.10.and.itype(i).ne.ntyp1) then
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.lt.3) then
do j=1,3
incr(j)=d_t(j,nres+i)
enddo
-!c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
- KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+&
+! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
incr(3)*incr(3))
endif
enddo
111 continue
!c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
!c & ' KEr_sc', KEr_sc
- KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
+ KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
!c write (iout,*) "KE_total",KE_total
return
end subroutine kinetic
implicit none
double precision KE_total
- integer i,j,k,iti,ichain,innt,inct
+ integer i,j,k,iti,ichain,innt,inct,mnum
double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
mag1,mag2,v(3)
#ifdef FIVEDIAG
!c & " inct",inct
do i=innt,inct-1
+ mnum=molnum(i)
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
!c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
do j=1,3
v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
enddo
!c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
- KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
enddo
!c write(iout,*) 'KEt_p', KEt_p
!c The translational part for the side chain virtual bond
!c Only now we can initialize incr with zeros. It must be equal
!c to the velocities of the first Calpha.
do i=innt,inct
- iti=iabs(itype(i))
- if (iti.eq.10) then
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then
!c write (iout,*) i,iti,(d_t(j,i),j=1,3)
do j=1,3
v(j)=d_t(j,i)
endif
!c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
!c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
- KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
enddo
!c goto 111
!c write(iout,*) 'KEt_sc', KEt_sc
!c The part due to stretching and rotation of the peptide groups
do i=innt,inct-1
+ mnum=molnum(i)
do j=1,3
incr(j)=d_t(j,i+1)-d_t(j,i)
enddo
+ if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
!c write (iout,*) i,(incr(j),j=1,3)
!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
- KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2)
- & +incr(3)*incr(3))
+ KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)&
+ +incr(3)*incr(3))
enddo
!c goto 111
!c write(iout,*) 'KEr_p', KEr_p
!c The rotational part of the side chain virtual bond
do i=innt,inct
- iti=iabs(itype(i))
- if (iti.ne.10) then
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+! if (iti.ne.10.and.mnum.lt.3) then
+ if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
do j=1,3
incr(j)=d_t(j,nres+i)-d_t(j,i)
enddo
!c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
- KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+
- & incr(3)*incr(3))
+ KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
+ incr(3)*incr(3))
endif
enddo
111 continue
!c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
!c & ' KEr_sc', KEr_sc
- KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
+ KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
!c write (iout,*) "KE_total",KE_tota
#else
write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
tt0=tcpu()
#endif
do itime=1,n_timestep
+ if (large) print *,itime,ntwe
if (ovrtim()) exit
if (large.and. mod(itime,ntwe).eq.0) &
write (iout,*) "itime",itime
enddo
#endif
else if (lang.eq.1 .or. lang.eq.4) then
+ print *,"before setup_fricmat"
call setup_fricmat
+ print *,"after setup_fricmat"
endif
write (iout,'(a,i10)') &
"Friction matrix reset based on surface area, itime",itime
call RESPA_step(itime)
else
! Variable time step algorithm.
+ if (large) print *,"before verlet_step"
call velverlet_step(itime)
+ if (large) print *,"after verlet_step"
endif
else
#ifdef BROWN
itime_mat=itime
if (ntwe.ne.0) then
if (mod(itime,ntwe).eq.0) then
- call returnbox
+! call returnbox
call statout(itime)
! call returnbox
! call check_ecartint
enddo
do i=nnt,nct
mnum=molnum(i)
- if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.lt.4) then
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then
do j=1,3
ind=ind+1
v_work(ind)=d_t(j,i+nres)
icount_scale=0
if (lang.eq.1) then
call sddir_precalc
+ if (large) print *,"after sddir_precalc"
else if (lang.eq.2 .or. lang.eq.3) then
#ifndef LANG0
call stochastic_force(stochforcvec)
call chainbuild_cart
if (rattle) call rattle1
if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
+ if (large) then !.and. mod(itime,ntwe).eq.0) then
write (iout,*) "Cartesian and internal coordinates: step 1"
call cartprint
call intout
!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
!el common /stochcalc/ stochforcvec
real(kind=8) :: time00
+ integer :: i
!
! Compute friction and stochastic forces
!
#ifdef MPI
time00=MPI_Wtime()
+ if (large) print *,"before friction_force"
call friction_force
+ if (large) print *,"after friction_force"
time_fric=time_fric+MPI_Wtime()-time00
time00=MPI_Wtime()
call stochastic_force(stochforcvec)
#ifdef FIVEDIAG
call fivediaginv_mult(dimen,fric_work, d_af_work)
call fivediaginv_mult(dimen,stochforcvec, d_as_work)
+ if (large) then
+ write(iout,*),"dimen",dimen
+ do i=1,dimen
+ write(iout,*) "fricwork",fric_work(i), d_af_work(i)
+ write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i)
+ enddo
+ endif
#else
call ginv_mult(fric_work, d_af_work)
call ginv_mult(stochforcvec, d_as_work)
do j=1,3
adt=(d_a_old(j,0)+d_af_work(j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
+! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time
dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
d_t(j,0)=d_t_old(j,0)+adt
do j=1,3
adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j)
dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
d_t(j,i)=d_t_old(j,i)+adt
do j=1,3
adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
-! write(iout,*) "adt",adt,"ads",adt2
+! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j)
dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
d_t(j,inres)=d_t_old(j,inres)+adt
ind=ind+3
endif
enddo
+
return
end subroutine sddir_verlet1
!-----------------------------------------------------------------------------
! forces (d_as_work)
!
#ifdef FIVEDIAG
- call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
+ call fivediaginv_mult(6*nres,stochforcvec, d_as_work1)
#else
call ginv_mult(stochforcvec, d_as_work1)
#endif
#endif
potE=potEcomp(0)
call cartgrad
+ write(iout,*) "before lagrangian"
call lagrangian
+ write(iout,*) "before max_accel"
call max_accel
if (amax*d_time .gt. dvmax) then
d_time=d_time*dvmax/amax
#ifdef FIVEDIAG
integer ichain,n,innt,inct,ibeg,ierr
double precision work(48*nres)
- integer iwork(maxres6)
- double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
- & Gvec(maxres2_chain,maxres2_chain)
- common /przechowalnia/Ghalf,Geigen,Gvec
+ integer iwork(6*nres)
+! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
+! Gvec(maxres2_chain,maxres2_chain)
+! common /przechowalnia/Ghalf,Geigen,Gvec
#ifdef DEBUG
double precision inertia(maxres2_chain,maxres2_chain)
#endif
#endif
-!#define DEBUG
+#define DEBUG
#ifdef FIVEDIAG
real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
- real(kind=8) :: sumx
+ real(kind=8) :: sumx,Ek2,Ek3,aux
#ifdef DEBUG
real(kind=8) ,allocatable, dimension(:) :: rsold
- real (kind=8),allocatable,dimension(:,:) :: matold
+ real (kind=8),allocatable,dimension(:,:) :: matold,inertia
integer :: iti
#endif
-#endif
integer :: i,j,ii,k,ind,mark,imark,mnum
! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
! First generate velocities in the eigenspace of the G matrix
! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
! call flush(iout)
- xv=0.0d0
- ii=0
- do i=1,dimen
- do k=1,3
- ii=ii+1
- sigv=dsqrt((Rb*t_bath)/geigen(i))
- lowb=-5*sigv
- highb=5*sigv
- d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+!#define DEBUG
#ifdef DEBUG
- write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
- " d_t_work_new",d_t_work_new(ii)
+ write (iout,*) "Random_vel, fivediag"
+ allocate(inertia(2*nres,2*nres))
#endif
- enddo
- enddo
+ d_t=0.0d0
+ Ek2=0.0d0
+ EK=0.0d0
+ Ek3=0.0d0
#ifdef DEBUG
-! diagnostics
- Ek1=0.0d0
- ii=0
- do i=1,dimen
- do k=1,3
- ii=ii+1
- Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
- enddo
- enddo
- write (iout,*) "Ek from eigenvectors",Ek1
- write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
-! end diagnostics
+ write(iout,*), nchain
#endif
-#ifdef FIVEDIAG
-! Transform velocities to UNRES coordinate space
- allocate (DL1(2*nres))
- allocate (DDU1(2*nres))
- allocate (DL2(2*nres))
- allocate (DDU2(2*nres))
- allocate (xsolv(2*nres))
- allocate (DML(2*nres))
- allocate (rs(2*nres))
+ do ichain=1,nchain
+ ind=0
+ ghalf=0.0d0
+ n=dimen_chain(ichain)
+ innt=iposd_chain(ichain)
+ inct=innt+n-1
#ifdef DEBUG
- allocate (rsold(2*nres))
- allocate (matold(2*nres,2*nres))
- matold=0
- matold(1,1)=DMorig(1)
- matold(1,2)=DU1orig(1)
- matold(1,3)=DU2orig(1)
- write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
-
- do i=2,dimen-1
- do j=1,dimen
- if (i.eq.j) then
- matold(i,j)=DMorig(i)
- matold(i,j-1)=DU1orig(i-1)
- if (j.ge.3) then
- matold(i,j-2)=DU2orig(i-2)
- endif
-
- endif
-
- enddo
- do j=1,dimen-1
- if (i.eq.j) then
- matold(i,j+1)=DU1orig(i)
-
- end if
- enddo
- do j=1,dimen-2
- if(i.eq.j) then
- matold(i,j+2)=DU2orig(i)
- endif
- enddo
- enddo
- matold(dimen,dimen)=DMorig(dimen)
- matold(dimen,dimen-1)=DU1orig(dimen-1)
- matold(dimen,dimen-2)=DU2orig(dimen-2)
- write(iout,*) "old gmatrix"
- call matout(dimen,dimen,2*nres,2*nres,matold)
+ write (iout,*) "Chain",ichain," n",n," start",innt
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),&
+ DU2orig(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
+ endif
+ enddo
#endif
- d_t_work=0.0d0
- do i=1,dimen
-! Find the ith eigenvector of the pentadiagonal inertiq matrix
- do imark=dimen,1,-1
- do j=1,imark-1
- DML(j)=DMorig(j)-geigen(i)
- enddo
- do j=imark+1,dimen
- DML(j-1)=DMorig(j)-geigen(i)
- enddo
- do j=1,imark-2
- DDU1(j)=DU1orig(j)
- enddo
- DDU1(imark-1)=DU2orig(imark-1)
- do j=imark+1,dimen-1
- DDU1(j-1)=DU1orig(j)
- enddo
- do j=1,imark-3
- DDU2(j)=DU2orig(j)
- enddo
- DDU2(imark-2)=0.0d0
- DDU2(imark-1)=0.0d0
- do j=imark,dimen-3
- DDU2(j)=DU2orig(j+1)
- enddo
- do j=1,dimen-3
- DL2(j+2)=DDU2(j)
- enddo
- do j=1,dimen-2
- DL1(j+1)=DDU1(j)
- enddo
-#ifdef DEBUG
- write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
- write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
- write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
- write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
- write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
- write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
-#endif
- rs=0.0d0
- if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
- if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
- if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
- if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
-#ifdef DEBUG
- rsold=rs
-#endif
-! write (iout,*) "Vector rs"
-! do j=1,dimen-1
-! write (iout,*) j,rs(j)
-! enddo
- xsolv=0.0d0
- call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
-
- if (mark.eq.1) then
-
-#ifdef DEBUG
-! Check solution
- do j=1,imark-1
- sumx=-geigen(i)*xsolv(j)
- do k=1,imark-1
- sumx=sumx+matold(j,k)*xsolv(k)
- enddo
- do k=imark+1,dimen
- sumx=sumx+matold(j,k)*xsolv(k-1)
- enddo
- write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
- enddo
- do j=imark+1,dimen
- sumx=-geigen(i)*xsolv(j-1)
- do k=1,imark-1
- sumx=sumx+matold(j,k)*xsolv(k)
- enddo
- do k=imark+1,dimen
- sumx=sumx+matold(j,k)*xsolv(k-1)
- enddo
- write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
- enddo
-! end check
- write (iout,*)&
- "Solution of equations system",i," for eigenvalue",geigen(i)
- do j=1,dimen-1
- write(iout,'(i5,f10.5)') j,xsolv(j)
- enddo
-#endif
- do j=dimen-1,imark,-1
- xsolv(j+1)=xsolv(j)
- enddo
- xsolv(imark)=1.0d0
-#ifdef DEBUG
- write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
- do j=1,dimen
- write(iout,'(i5,f10.5)') j,xsolv(j)
- enddo
-#endif
-! Normalize ith eigenvector
- sumx=0.0d0
- do j=1,dimen
- sumx=sumx+xsolv(j)**2
- enddo
- sumx=dsqrt(sumx)
- do j=1,dimen
- xsolv(j)=xsolv(j)/sumx
- enddo
+ ghalf(ind+1)=dmorig(innt)
+ ghalf(ind+2)=du1orig(innt)
+ ghalf(ind+3)=dmorig(innt+1)
+ ind=ind+3
+ do i=3,n
+ ind=ind+i-3
+ write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
+ " indu1",innt+i-1," indm",innt+i
+ ghalf(ind+1)=du2orig(innt-1+i-2)
+ ghalf(ind+2)=du1orig(innt-1+i-1)
+ ghalf(ind+3)=dmorig(innt-1+i)
+!c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
+!c & "DU1",innt-1+i-1,"DM ",innt-1+i
+ ind=ind+3
+ enddo
#ifdef DEBUG
- write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
- do j=1,dimen
- write(iout,'(i5,3f10.5)') j,xsolv(j)
- enddo
+ ind=0
+ do i=1,n
+ do j=1,i
+ ind=ind+1
+ inertia(i,j)=ghalf(ind)
+ inertia(j,i)=ghalf(ind)
+ enddo
+ enddo
#endif
-! All done at this point for eigenvector i, exit loop
- exit
-
- endif ! mark.eq.1
-
- enddo ! imark
-
- if (mark.ne.1) then
- write (iout,*) "Unable to find eigenvector",i
- endif
-
-! write (iout,*) "i=",i
- do k=1,3
-! write (iout,*) "k=",k
- do j=1,dimen
- ind=(j-1)*3+k
-! write(iout,*) "ind",ind," ind1",3*(i-1)+k
- d_t_work(ind)=d_t_work(ind) &
- +xsolv(j)*d_t_work_new(3*(i-1)+k)
- enddo
- enddo
- enddo ! i (loop over the eigenvectors)
-
#ifdef DEBUG
- write (iout,*) "d_t_work"
- do i=1,3*dimen
- write (iout,"(i5,f10.5)") i,d_t_work(i)
- enddo
- Ek1=0.0d0
- ii=0
- do i=nnt,nct
-! if (itype(i,1).eq.10) then
- mnum=molnum(i)
- if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
- iti=iabs(itype(i,mnum))
-! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
- if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
- .or.(mnum.ge.4)) then
- j=ii+3
- else
- j=ii+6
+ write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
+ write (iout,*) "Five-diagonal inertia matrix, lower triangle"
+! call matoutr(n,ghalf)
+#endif
+ call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
+ if (large) then
+ write (iout,'(//a,i3)')&
+ "Eigenvectors and eigenvalues of the G matrix chain",ichain
+ call eigout(n,n,nres*2,nres*2,Gvec,Geigen)
endif
- if (i.lt.nct) then
+#ifdef DIAGCHECK
+!c check diagonalization
+ do i=1,n
+ do j=1,n
+ aux=0.0d0
+ do k=1,n
+ do l=1,n
+ aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
+ enddo
+ enddo
+ if (i.eq.j) then
+ write (iout,*) i,j,aux,geigen(i)
+ else
+ write (iout,*) i,j,aux
+ endif
+ enddo
+ enddo
+#endif
+ xv=0.0d0
+ ii=0
+ do i=1,n
do k=1,3
-! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
- Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
- 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
+ ii=ii+1
+ sigv=dsqrt((Rb*t_bath)/geigen(i))
+ lowb=-5*sigv
+ highb=5*sigv
+ d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+ EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
+!c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
+!c & " d_t_work_new",d_t_work_new(ii)
enddo
- endif
- if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.lt.4)) ii=ii+3
- write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
- write (iout,*) "ii",ii
+ enddo
do k=1,3
- ii=ii+1
- write (iout,*) "k",k," ii",ii,"EK1",EK1
- if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.lt.4))&
- Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
- Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
+ do i=1,n
+ ind=(i-1)*3+k
+ d_t_work(ind)=0.0d0
+ do j=1,n
+ d_t_work(ind)=d_t_work(ind)&
+ +Gvec(i,j)*d_t_work_new((j-1)*3+k)
+ enddo
+!c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+!c call flush(iout)
+ enddo
enddo
- write (iout,*) "i",i," ii",ii
- enddo
- write (iout,*) "Ek from d_t_work",Ek1
- write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
-#endif
- if(allocated(DDU1)) deallocate(DDU1)
- if(allocated(DDU2)) deallocate(DDU2)
- if(allocated(DL2)) deallocate(DL2)
- if(allocated(DL1)) deallocate(DL1)
- if(allocated(xsolv)) deallocate(xsolv)
- if(allocated(DML)) deallocate(DML)
- if(allocated(rs)) deallocate(rs)
#ifdef DEBUG
- if(allocated(matold)) deallocate(matold)
- if(allocated(rsold)) deallocate(rsold)
-#endif
- ind=1
- do j=nnt,nct
+ aux=0.0d0
do k=1,3
- d_t(k,j)=d_t_work(ind)
- ind=ind+1
+ do i=1,n
+ do j=1,n
+ aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
+ enddo
+ enddo
enddo
- mnum=molnum(j)
- if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.lt.4)) then
- do k=1,3
- d_t(k,j+nres)=d_t_work(ind)
+ Ek3=Ek3+aux/2
+#endif
+!c Transfer to the d_t vector
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ ind=0
+!c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+ do i=innt,inct
+ do j=1,3
ind=ind+1
+ d_t(j,i)=d_t_work(ind)
enddo
- end if
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ do j=1,3
+ ind=ind+1
+ d_t(j,i+nres)=d_t_work(ind)
+ enddo
+ endif
+ enddo
enddo
-#ifdef DEBUG
- write (iout,*) "Random velocities in the Calpha,SC space"
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
+ if (large) then
+ write (iout,*)
+ write (iout,*) "Random velocities in the Calpha,SC space"
+ do i=1,nres
+ mnum=molnum(i)
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
+ restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ call kinetic_CASC(Ek1)
+!
+! Transform the velocities to virtual-bond space
+!
+#define WLOS
+#ifdef WLOS
+ if (nnt.eq.1) then
+ d_t(:,0)=d_t(:,1)
+ endif
+ do i=1,nres
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then
+ do j=1,3
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ else
+ do j=1,3
+ d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ end if
enddo
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+ d_t(:,nres)=0.0d0
+ d_t(:,nct)=0.0d0
+ d_t(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ d_t(:,0)=d_t(:,1)
+ d_t(:,1)=0.0d0
+ endif
+!c d_a(:,0)=d_a(:,1)
+!c d_a(:,1)=0.0d0
+!c write (iout,*) "Shifting accelerations"
+ do ichain=2,nchain
+!c write (iout,*) "ichain",chain_border1(1,ichain)-1,
+!c & chain_border1(1,ichain)
+ d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
+ d_t(:,chain_border1(1,ichain))=0.0d0
+ enddo
+!c write (iout,*) "Adding accelerations"
+ do ichain=2,nchain
+!c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+!c & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=&
+ d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+ do ichain=2,nchain
+ write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
+ chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=&
+ d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
enddo
-#endif
+#else
+ ibeg=0
+!c do j=1,3
+!c d_t(j,0)=d_t(j,nnt)
+!c enddo
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+!c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+!c write (iout,*) "ibeg",ibeg
do j=1,3
- d_t(j,0)=d_t(j,nnt)
+ d_t(j,ibeg)=d_t(j,innt)
enddo
- do i=nnt,nct
-! if (itype(i,1).eq.10) then
- mnum=molnum(i)
- if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
- .or.(mnum.ge.4)) then
+ ibeg=inct+1
+ do i=innt,inct
+ mnum=molnum(i)
+ if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
+!c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
do j=1,3
d_t(j,i)=d_t(j,i+1)-d_t(j,i)
enddo
enddo
end if
enddo
+ enddo
+#endif
+ if (large) then
+ write (iout,*)
+ write (iout,*)&
+ "Random velocities in the virtual-bond-vector space"
+ write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+ do i=1,nres
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
+ restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ write (iout,*)
+ write (iout,*) "Kinetic energy from inertia matrix eigenvalues",&
+ Ek
+ write (iout,*)&
+ "Kinetic temperatures from inertia matrix eigenvalues",&
+ 2*Ek/(3*dimen*Rb)
#ifdef DEBUG
- write (iout,*)"Random velocities in the virtual-bond-vector space"
- do i=nnt,nct-1
- write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
+ write (iout,*) "Kinetic energy from inertia matrix",Ek3
+ write (iout,*) "Kinetic temperatures from inertia",&
+ 2*Ek3/(3*dimen*Rb)
+#endif
+ write (iout,*) "Kinetic energy from velocities in CA-SC space",&
+ Ek1
+ write (iout,*)&
+ "Kinetic temperatures from velovities in CA-SC space",&
+ 2*Ek1/(3*dimen*Rb)
+ call kinetic(Ek1)
+ write (iout,*)&
+ "Kinetic energy from virtual-bond-vector velocities",Ek1
+ write (iout,*)&
+ "Kinetic temperature from virtual-bond-vector velocities ",&
+ 2*Ek1/(dimen3*Rb)
+ endif
+#else
+ xv=0.0d0
+ ii=0
+ do i=1,dimen
+ do k=1,3
+ ii=ii+1
+ sigv=dsqrt((Rb*t_bath)/geigen(i))
+ lowb=-5*sigv
+ highb=5*sigv
+ d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+#ifdef DEBUG
+ write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
+ " d_t_work_new",d_t_work_new(ii)
+#endif
+ enddo
enddo
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+#ifdef DEBUG
+! diagnostics
+ Ek1=0.0d0
+ ii=0
+ do i=1,dimen
+ do k=1,3
+ ii=ii+1
+ Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
+ enddo
enddo
- call kinetic(Ek1)
- write (iout,*) "Ek from d_t_work",Ek1
+ write (iout,*) "Ek from eigenvectors",Ek1
write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
+! end diagnostics
#endif
-#else
+
do k=0,2
do i=1,dimen
ind=(i-1)*3+k+1
! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
! call flush(iout)
! write(iout,*) "end init MD"
+#undef DEBUG
return
end subroutine random_vel
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! moments.f
!-----------------------------------------------------------------------------
+#ifdef FIVEDIAG
+ subroutine inertia_tensor
+ use comm_gucio
+ use energy_data
+ real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,&
+ eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),&
+ vpp(3,0:MAXRES),vs_p(3),pr1(3,3),&
+ pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP
+ integer iti,inres,i,j,k,mnum,mnum1
+ do i=1,3
+ do j=1,3
+ Im(i,j)=0.0d0
+ pr1(i,j)=0.0d0
+ pr2(i,j)=0.0d0
+ enddo
+ L(i)=0.0d0
+ cm(i)=0.0d0
+ vrot(i)=0.0d0
+ enddo
+ M_PEP=0.0d0
+
+!c caulating the center of the mass of the protein
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ M_PEP=M_PEP+mp(mnum)
+
+ do j=1,3
+ cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
+ enddo
+ enddo
+! do j=1,3
+! cm(j)=mp*cm(j)
+! enddo
+ M_SC=0.0d0
+ do i=nnt,nct
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ iti=iabs(itype(i,mnum))
+ if (iti.eq.ntyp1_molec(mnum)) cycle
+ M_SC=M_SC+msc(iabs(iti),mnum)
+ inres=i+nres
+ do j=1,3
+ cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/(M_SC+M_PEP)
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (iti.eq.ntyp1_molec(mnum)) cycle
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ iti=iabs(itype(i,mnum))
+ if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ inres=i+nres
+ Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
+ dc_norm(1,inres))*vbld(inres)*vbld(inres)
+ Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ endif
+ enddo
+
+ call angmom(cm,L)
+ Im(2,1)=Im(1,2)
+ Im(3,1)=Im(1,3)
+ Im(3,2)=Im(2,3)
+
+!c Copng the Im matrix for the djacob subroutine
+ do i=1,3
+ do j=1,3
+ Imcp(i,j)=Im(i,j)
+ Id(i,j)=0.0d0
+ enddo
+ enddo
+!c Finding the eigenvectors and eignvalues of the inertia tensor
+ call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
+ do i=1,3
+ if (dabs(eigval(i)).gt.1.0d-15) then
+ Id(i,i)=1.0d0/eigval(i)
+ else
+ Id(i,i)=0.0d0
+ endif
+ enddo
+ do i=1,3
+ do j=1,3
+ Imcp(i,j)=eigvec(j,i)
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ do k=1,3
+ pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
+ enddo
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ do k=1,3
+ pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
+ enddo
+ enddo
+ enddo
+!c Calculating the total rotational velocity of the molecule
+ do i=1,3
+ do j=1,3
+ vrot(i)=vrot(i)+pr2(i,j)*L(j)
+ enddo
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
+ if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
+ .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
+ itype(i,mnum).ne.ntyp1_molec(mnum) &
+ .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ call vecpr(vrot(1),dc(1,i),vp)
+ do j=1,3
+ d_t(j,i)=d_t(j,i)-vp(j)
+ enddo
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
+ inres=i+nres
+ call vecpr(vrot(1),dc(1,inres),vp)
+ do j=1,3
+ d_t(j,inres)=d_t(j,inres)-vp(j)
+ enddo
+ endif
+ enddo
+ call angmom(cm,L)
+ return
+ end subroutine inertia_tensor
+!------------------------------------------------------------
+ subroutine angmom(cm,L)
+ implicit none
+ double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
+ pp(3),mscab
+ integer iti,inres,i,j,mnum,mnum1
+!c Calculate the angular momentum
+ do j=1,3
+ L(j)=0.0d0
+ enddo
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ call vecpr(pr(1),v(1),vp)
+ do j=1,3
+ L(j)=L(j)+mp(mnum)*vp(j)
+ enddo
+ do j=1,3
+ pr(j)=0.5d0*dc(j,i)
+ pp(j)=0.5d0*d_t(j,i)
+ enddo
+ call vecpr(pr(1),pp(1),vp)
+ do j=1,3
+ L(j)=L(j)+Ip(mnum)*vp(j)
+ enddo
+ enddo
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
+ do j=1,3
+ v(j)=incr(j)+d_t(j,inres)
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)
+ enddo
+ endif
+ call vecpr(pr(1),v(1),vp)
+!c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
+!c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
+ if (mnum.gt.4) then
+ mscab=0.0d0
+ else
+ mscab=msc(iti,mnum)
+ endif
+ do j=1,3
+ L(j)=L(j)+mscab*vp(j)
+ enddo
+!c write (iout,*) "L",(l(j),j=1,3)
+ if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
+ do j=1,3
+ v(j)=incr(j)+d_t(j,inres)
+ enddo
+ call vecpr(dc(1,inres),d_t(1,inres),vp)
+ do j=1,3
+ L(j)=L(j)+Isc(iti,mnum)*vp(j)
+ enddo
+ endif
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+ return
+ end subroutine angmom
+!---------------------------------------------------------------
+ subroutine vcm_vel(vcm)
+ double precision vcm(3),vv(3),summas,amas
+ integer i,j,mnum,mnum1
+ do j=1,3
+ vcm(j)=0.0d0
+ vv(j)=d_t(j,0)
+ enddo
+ summas=0.0d0
+ do i=nnt,nct
+ mnum=molnum(i)
+ if ((mnum.ge.5).or.(mnum.eq.3))&
+ mp(mnum)=msc(itype(i,mnum),mnum)
+ if (i.lt.nct) then
+ summas=summas+mp(mnum)
+ do j=1,3
+ vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
+ enddo
+ endif
+ if (mnum.ne.4) then
+ amas=msc(iabs(itype(i,mnum)),mnum)
+ else
+ amas=0.0d0
+ endif
+! amas=msc(iabs(itype(i)))
+ summas=summas+amas
+! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.lt.4)) then
+ do j=1,3
+ vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
+ enddo
+ else
+ do j=1,3
+ vcm(j)=vcm(j)+amas*vv(j)
+ enddo
+ endif
+ do j=1,3
+ vv(j)=vv(j)+d_t(j,i)
+ enddo
+ enddo
+!c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+ do j=1,3
+ vcm(j)=vcm(j)/summas
+ enddo
+ return
+ end subroutine vcm_vel
+#else
subroutine inertia_tensor
! Calculating the intertia tensor for the entire protein in order to
enddo
return
end subroutine vcm_vel
+#endif
!-----------------------------------------------------------------------------
! rattle.F
!-----------------------------------------------------------------------------
!el common /syfek/ gamvec
#ifdef FIVEDIAG
integer iposc,ichain,n,innt,inct
- double precision rs(nres+2)
+ double precision rs(nres*2)
+ real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
+#else
+ real(kind=8) :: v_work(6*nres)
#endif
- real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
+ real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,&
!el ginvfric(2*nres,2*nres) !maxres2=2*maxres
!el common /przechowalnia/ ginvfric
! if (large) lprn=.true.
! if (large) checkmode=.true.
#ifdef FIVEDIAG
-c Here accelerations due to friction forces are computed right after forces.
+!c Here accelerations due to friction forces are computed right after forces.
d_t_work(:6*nres)=0.0d0
do j=1,3
v_work(j,1)=d_t(j,0)
do ichain=1,nchain
n=dimen_chain(ichain)
iposc=iposd_chain(ichain)
-c write (iout,*) "friction_force j",j," ichain",ichain,
-c & " n",n," iposc",iposc,iposc+n-1
+!c write (iout,*) "friction_force j",j," ichain",ichain,
+!c & " n",n," iposc",iposc,iposc+n-1
innt=chain_border(1,ichain)
inct=chain_border(2,ichain)
-c diagnostics
-c innt=chain_border(1,1)
-c inct=chain_border(2,1)
+!c diagnostics
+!c innt=chain_border(1,1)
+!c inct=chain_border(2,1)
do i=innt,inct
vvec(ind+1)=v_work(j,i)
ind=ind+1
write (iout,*) "vvec ind",ind," n",n
write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
#endif
-c write (iout,*) "chain",i," ind",ind," n",n
+!c write (iout,*) "chain",i," ind",ind," n",n
#ifdef TIMING
#ifdef MPI
time01=MPI_Wtime()
time01=tcpu()
#endif
#endif
- call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
- & DU2fric(iposc),vvec(iposc),rs)
+! if (large) print *,"before fivediagmult"
+ call fivediagmult(n,DMfric(iposc),DU1fric(iposc),&
+ DU2fric(iposc),vvec(iposc),rs)
+! if (large) print *,"after fivediagmult"
+
#ifdef TIMING
#ifdef MPI
time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
write (iout,'(f10.5)') (rs(i),i=1,n)
#endif
do i=iposc,iposc+n-1
-c write (iout,*) "ichain",ichain," i",i," j",j,
-c & "index",3*(i-1)+j,"rs",rs(i-iposc+1)
+! write (iout,*) "ichain",ichain," i",i," j",j,&
+! "index",3*(i-1)+j,"rs",rs(i-iposc+1)
fric_work(3*(i-1)+j)=-rs(i-iposc+1)
enddo
enddo
!#endif
! include 'COMMON.IOUNITS'
integer :: IERROR
- integer :: i,j,ind,ind1,m
+ integer :: i,j,ind,ind1,m,ichain,innt,inct
logical :: lprn = .false.
real(kind=8) :: dtdi !el ,gamvec(2*nres)
!el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
nres2=2*nres
nres6=6*nres
#ifdef MPI
+#ifndef FIVEDIAG
if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
if (fg_rank.ne.king) goto 10
#endif
+#endif
! nres2=2*nres
! nres6=6*nres
if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
+#ifndef FIVEDIAG
if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
!el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
+#endif
+#ifdef FIVEDIAG
+ if (.not.allocated(DMfric)) allocate(DMfric(nres2))
+ if (.not.allocated(DU1fric)) allocate(DU1fric(nres2))
+ if (.not.allocated(DU2fric)) allocate(DU2fric(nres2))
+! Load the friction coefficients corresponding to peptide groups
+ ind1=0
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ ind1=ind1+1
+ gamvec(ind1)=gamp(mnum)
+ enddo
+!HERE TEST
+ if (molnum(nct).eq.5) then
+ mnum=molnum(i)
+ ind1=ind1+1
+ gamvec(ind1)=gamp(mnum)
+ endif
+! Load the friction coefficients corresponding to side chains
+ m=nct-nnt
+ ind=0
+ do j=1,2
+ gamsc(ntyp1_molec(j),j)=1.0d0
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ ind=ind+1
+ ii = ind+m
+ iti=itype(i,mnum)
+ gamvec(ii)=gamsc(iabs(iti),mnum)
+ enddo
+ if (surfarea) call sdarea(gamvec)
+ DMfric=0.0d0
+ DU1fric=0.0d0
+ DU2fric=0.0d0
+ ind=1
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+!c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+!c DMfric part
+ mnum=molnum(innt)
+ DMfric(ind)=gamvec(innt-nnt+1)/4
+ if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then
+ DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(m+innt-nnt+1)
+ ind=ind+2
+ endif
+!c write (iout,*) "DMfric init ind",ind
+!c DMfric
+ do i=innt+1,inct-1
+ mnum=molnum(i)
+ DMfric(ind)=gamvec(i-nnt+1)/2
+ if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then
+ DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(m+i-nnt+1)
+ ind=ind+2
+ endif
+ enddo
+!c write (iout,*) "DMfric endloop ind",ind
+ if (inct.gt.innt) then
+ DMfric(ind)=gamvec(inct-1-nnt+1)/4
+ mnum=molnum(inct)
+ if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then
+ DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(inct+m-nnt+1)
+ ind=ind+2
+ endif
+ endif
+!c write (iout,*) "DMfric end ind",ind
+ enddo
+!c DU1fric part
+ do ichain=1,nchain
+ mnum=molnum(i)
+
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct
+ if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
+ ind=ind+2
+ else
+ DU1fric(ind)=gamvec(i-nnt+1)/4
+ ind=ind+1
+ endif
+ enddo
+ enddo
+!c DU2fric part
+ do ichain=1,nchain
+ mnum=molnum(i)
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct-1
+ if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
+ DU2fric(ind)=gamvec(i-nnt+1)/4
+ DU2fric(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU2fric(ind)=0.0d0
+ ind=ind+1
+ endif
+ enddo
+ enddo
+ if (lprn) then
+ write(iout,*)"The upper part of the five-diagonal friction matrix"
+ do ichain=1,nchain
+ write (iout,'(a,i5)') 'Chain',ichain
+ innt=iposd_chain(ichain)
+ inct=iposd_chain(ichain)+dimen_chain(ichain)-1
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),&
+ DU2fric(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
+ endif
+ enddo
+ enddo
+ endif
+ 10 continue
+#else
+
+
! Zeroing out fricmat
do i=1,dimen
do j=1,dimen
! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
endif
#endif
+#endif
return
end subroutine setup_fricmat
!-----------------------------------------------------------------------------
do j=1,3
stochforcvec(ind+j)=0.5d0*force(j,innt)
enddo
- if (iabs(itype(innt),molnum(iint)).eq.10) then
+ if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
do j=1,3
stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
enddo
do j=1,3
stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
enddo
- if (iabs(itype(i,molnum(i)).eq.10) then
+ if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then
do j=1,3
stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
enddo
do j=1,3
stochforcvec(ind+j)=0.5d0*force(j,inct-1)
enddo
- if (iabs(itype(inct),molnum(inct)).eq.10) then
+ if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then
do j=1,3
stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
enddo
#ifdef FIVEDIAG
allocate(DM(nres2),DU1(nres2),DU2(nres2))
allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
+ allocate(Gvec(nres2,nres2))
#else
write (iout,*) "Before A Allocation",nfgtasks-1
call flush(iout)
! include 'DIMENSIONS'
use comm_cipiszcze
use energy_data
+ use energy, only: grad_transform
use geometry_data, only: nres
use control_data !el, only: mucadyn,lmuca
#ifdef MPI
! include 'COMMON.CONTROL'
! include 'COMMON.MUCA'
! include 'COMMON.TIME1'
- integer ::i,j,ind,itime,mnum
+ integer ::i,j,ind,itime,mnum,innt,inct,inct_prev,ichain,n,mark
real(kind=8) :: zapas(6*nres) !,muca_factor !maxres6=6*maxres
- real(kind=8) :: rs(dimen),xsolv(dimen)
+! the line below might be wrong
+#ifdef FIVEDIAG
+ real(kind=8) :: rs(2*nres),xsolv(2*nres)
#ifdef CHECK5DSOL
- real(kind=8) :: rscheck(dimen),rsold(dimen)
+ real(kind=8) :: rscheck(2*nres),rsold(2*nres)
#endif
- logical :: lprn = .false.
+#endif
+ logical :: lprn = .true.
!el common /cipiszcze/ itime
itime = itt_comm
time00=MPI_Wtime()
#endif
#ifdef FIVEDIAG
+ call grad_transform
+ d_a=0.0d0
if (lprn) then
write (iout,*) "Potential forces backbone"
- do i=nnt,nct
+ do i=1,nres
write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
enddo
write (iout,*) "Potential forces sidechain"
do i=nnt,nct
- mnum=molnum(i)
- if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.lt.4)&
- write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
+! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+ write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
enddo
endif
- do j=1,3
- ind=1
- do i=nnt,nct
- mnum=molnum(i)
- if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.ge.4)then
- rs(ind)=-gcart(j,i)-gxcart(j,i)
- ind=ind+1
- else
- rs(ind)=-gcart(j,i)
- rs(ind+1)=-gxcart(j,i)
- ind=ind+2
- end if
- enddo
+ do ichain=1,nchain
+ n=dimen_chain(ichain)
+ innt=iposd_chain(ichain)
+ do j=1,3
+ ind=1
+ do i=chain_border(1,ichain),chain_border(2,ichain)
+ mnum=molnum(i)
+ if (itype(i,1).eq.10.or.mnum.ge.3)then
+ rs(ind)=-gcart(j,i)-gxcart(j,i)
+ ind=ind+1
+ else
+ rs(ind)=-gcart(j,i)
+ rs(ind+1)=-gxcart(j,i)
+ ind=ind+2
+ end if
+ enddo
+
#ifdef CHECK5DSOL
rsold=rs
#endif
if (lprn) then
- write(iout,*) "RHS of the 5-diag equations system"
- do i=1,dimen
+ write(iout,*) "RHS of the 5-diag equations system",&
+ ichain," j",j
+ do i=1,n
write(iout,*) i,rs(i)
enddo
endif
- call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
+ call FDISYS (n,DM(innt),DU1(innt),DU2(innt),rs,xsolv)
if (lprn) then
write (iout,*) "Solution of the 5-diagonal equations system"
- do i=1,dimen
+ do i=1,n
write (iout,'(i5,f10.5)') i,xsolv(i)
enddo
endif
#ifdef CHECK5DSOL
! Check the solution
- rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
- DU2orig(1)*xsolv(3)
- rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
- DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
-
- do i=3,dimen-2
- rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
- xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
- xsolv(i+1)+DU2orig(i)*xsolv(i+2)
- enddo
- rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
- DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
- xsolv(dimen-1)+DU1orig(dimen-1)*&
- xsolv(dimen)
- rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
- xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
-
- do i=1,dimen
- write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
- "ratio",rscheck(i)/rsold(i)
- enddo
+ call fivediagmult(n,DMorig(innt),DU1orig(innt),DU2orig(innt),
+ & xsolv,rscheck)
+ do i=1,n
+ write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),
+ & "ratio",rscheck(i)/rsold(i)
+ enddo
! end check
#endif
- ind=1
- do i=nnt,nct
- mnum=molnum(i)
- if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or. mnum.ge.4)then
- d_a(j,i)=xsolv(ind)
- ind=ind+1
- else
- d_a(j,i)=xsolv(ind)
- d_a(j,i+nres)=xsolv(ind+1)
- ind=ind+2
- end if
- enddo
- enddo
- if (lprn) then
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
- enddo
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
- enddo
- endif
- do j=1,3
- d_a(j,0)=d_a(j,nnt)
- enddo
- do i=nnt,nct
- mnum=molnum(i)
-! if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
- if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or. mnum.ge.4)then
- do j=1,3
- d_a(j,i)=d_a(j,i+1)-d_a(j,i)
- enddo
- else
- do j=1,3
- d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
- d_a(j,i)=d_a(j,i+1)-d_a(j,i)
- enddo
-
-
- end if
- enddo
-
+#undef CHECK5DSOL
+ ind=1
+ do i=chain_border(1,ichain),chain_border(2,ichain)
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or.mnum.ge.3) then
+ d_a(j,i)=xsolv(ind)
+ ind=ind+1
+ else
+ d_a(j,i)=xsolv(ind)
+ d_a(j,i+nres)=xsolv(ind+1)
+ ind=ind+2
+ end if
+ enddo
+ enddo ! j
+ enddo ! ichain
if (lprn) then
- write(iout,*) 'acceleration 3D FIVEDIAG'
- write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
- do i=nnt,nct-1
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
+ write (iout,*) "Acceleration in CA and SC oordinates"
+ do i=1,nres
+ write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
enddo
do i=nnt,nct
- write (iout,'(i3,3f10.5,3x,3f10.5)') &
- i+nres,(d_a(j,i+nres),j=1,3)
+ write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+!C Conevert d_a to virtual-bon-vector basis
+#define WLOS
+#ifdef WLOS
+!c write (iout,*) "WLOS"
+ if (nnt.eq.1) then
+ d_a(:,0)=d_a(:,1)
+ endif
+ do i=1,nres
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum) .or.mnum.ge.3) then
+ do j=1,3
+ d_a(j,i)=d_a(j,i+1)-d_a(j,i)
+ enddo
+ else
+ do j=1,3
+ d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
+ d_a(j,i)=d_a(j,i+1)-d_a(j,i)
+ enddo
+ end if
+ enddo
+ d_a(:,nres)=0.0d0
+ d_a(:,nct)=0.0d0
+ d_a(:,2*nres)=0.0d0
+!c d_a(:,0)=d_a(:,1)
+!c d_a(:,1)=0.0d0
+!c write (iout,*) "Shifting accelerations"
+ if (nnt.gt.1) then
+ d_a(:,0)=d_a(:,1)
+ d_a(:,1)=0.0d0
+ endif
+#define CHUJ
+#ifdef CHUJ
+ do ichain=2,nchain
+!c write (iout,*) "ichain",chain_border1(1,ichain)-1,
+!c & chain_border1(1,ichain)
+ d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain))
+ d_a(:,chain_border1(1,ichain))=0.0d0
+ enddo
+!c write (iout,*) "Adding accelerations"
+ do ichain=2,nchain
+!c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+!c & chain_border(2,ichain-1)
+ d_a(:,chain_border1(1,ichain)-1)=&
+ d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
+ d_a(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+#endif
+#else
+ inct_prev=0
+ do j=1,3
+ aaux(j)=0.0d0
+ enddo
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do j=1,3
+ d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
+ enddo
+ inct_prev=inct+1
+ do i=innt,inct
+ if (itype(i).ne.10) then
+ do j=1,3
+ d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
+ enddo
+ endif
+ enddo
+ do j=1,3
+ aaux(j)=d_a(j,inct)
+ enddo
+ do i=innt,inct
+ do j=1,3
+ d_a(j,i)=d_a(j,i+1)-d_a(j,i)
+ enddo
+ enddo
+ enddo
+#endif
+ if (lprn) then
+ write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
+ enddo
+ do i=nnt,nct
+ write (iout,'(i3,3f10.5,3x,3f10.5)')&
+ i,(d_a(j,i+nres),j=1,3)
enddo
endif
+
#else
! Old procedure
do j=1,3
real(kind=8) :: coeff,mscab
integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
real(kind=8) :: ip4
- integer :: iz,mnum
+ integer :: iz,mnum,ichain,n,dimenp,innt,inct
print *,"just entered"
relfeh=1.0d-14
nres2=2*nres
print *,"FIVEDIAG"
write (iout,*) "before FIVEDIAG"
+ if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
#ifndef FIVEDIAG
write (iout,*) "ALLOCATE"
print *,"ALLOCATE"
if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
- if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
+! if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
if(.not.allocated(Bmat)) allocate(Bmat(nres2,nres2))
if(.not.allocated(matmult)) allocate(matmult(nres2,nres2))
#endif
!
! Determine the number of degrees of freedom (dimen) and the number of
! sites (dimen1)
- dimen=(nct-nnt+1)+nside
- dimen1=(nct-nnt)+(nct-nnt+1)
- dimen3=dimen*3
- write (iout,*) "nnt",nnt," nct",nct," nside",nside
+! dimen=(nct-nnt+1)+nside
+! dimen1=(nct-nnt)+(nct-nnt+1)
+! dimen3=dimen*3
+! write (iout,*) "nnt",nnt," nct",nct," nside",nside
#ifdef FIVEDIAG
#ifdef CRYST_BOND
ip4=ip/4
endif
enddo
#else
- ip4=ip(1)/4
- DM(1)=mp(1)/4+ip4
- mnum=molnum(nnt)
- if (iabs(itype(nnt,1)).eq.10) then
- DM(1)=DM(1)+msc(10,1)
- ind=2
- nind=1
- else
- DM(1)=DM(1)+isc(iabs(itype(nnt,mnum)),mnum)
- DM(2)=msc(iabs(itype(nnt,mnum)),mnum)+isc(iabs(itype(nnt,mnum)),mnum)
- ind=3
- nind=2
- endif
- do i=nnt+1,nct-1
-! if (iabs(itype(i,1)).eq.ntyp1) cycle
+ dimen=0
+ dimen1=0
+ dimenp=0
+ do ichain=1,nchain
+ dimen=dimen+chain_length(ichain)
+ dimen1=dimen1+2*chain_length(ichain)-1
+ dimenp=dimenp+chain_length(ichain)-1
+ enddo
+ write (iout,*) "Number of Calphas",dimen
+ write (iout,*) "Number of sidechains",nside
+ write (iout,*) "Number of peptide groups",dimenp
+ dimen=dimen+nside ! number of centers
+ dimen3=3*dimen ! degrees of freedom
+ write (iout,*) "Number of centers",dimen
+ write (iout,*) "Degrees of freedom:",dimen3
+! ip4=ip/4
+ ind=1
+ do ichain=1,nchain
+ iposd_chain(ichain)=ind
+ innt=chain_border(1,ichain)
+ mnum=molnum(innt)
+ inct=chain_border(2,ichain)
+ DM(ind)=mp(mnum)/4+ip(mnum)/4
+ if (iabs(itype(innt,1)).eq.10.or.molnum(innt).gt.2) then
+ DM(ind)=DM(ind)+msc(itype(innt,mnum),mnum)
+ ind=ind+1
+ nind=1
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(innt,mnum)),mnum)
+ DM(ind+1)=msc(iabs(itype(innt,mnum)),mnum)+isc(iabs(itype(innt,mnum)),mnum)
+ ind=ind+2
+ nind=2
+ endif
+ write (iout,*) "ind",ind," nind",nind
+ do i=innt+1,inct-1
mnum=molnum(i)
- DM(ind)=2*ip4+mp(1)/2
- if (iabs(itype(i,1)).eq.10 .or. &
- iabs(itype(i,mnum)).eq.ntyp1_molec(mnum) .or. mnum.ge.4) then
- if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10,1)
- if (mnum.eq.5) DM(ind)=DM(ind)+msc(itype(i,mnum),mnum)
- ind=ind+1
- else
- DM(ind)=DM(ind)+isc(iabs(itype(i,mnum)),mnum)
- DM(ind+1)=msc(iabs(itype(i,mnum)),mnum)+isc(iabs(itype(i,mnum)),mnum)
- ind=ind+2
- endif
- enddo
- if (nct.gt.nnt) then
- DM(ind)=ip4+mp(1)/4
- if (iabs(itype(nct,1)).eq.10) then
- DM(ind)=DM(ind)+msc(10,1)
- nind=ind
- else
- mnum=molnum(nct)
- DM(ind)=DM(ind)+isc(iabs(itype(nct,mnum)),mnum)
- DM(ind+1)=msc(iabs(itype(nct,mnum)),mnum)+isc(iabs(itype(nct,mnum)),mnum)
- nind=ind+1
- endif
+! if (iabs(itype(i)).eq.ntyp1) cycle
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
+ DM(ind)=2*ip(mnum)/4+mp(mnum)/2
+ if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2) then
+ if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2)&
+ DM(ind)=DM(ind)+msc(itype(i,molnum(i)),mnum)
+ ind=ind+1
+ nind=nind+1
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(i,mnum)),mnum)
+ DM(ind+1)=msc(iabs(itype(i,mnum)),mnum)+isc(iabs(itype(i,mnum)),mnum)
+ ind=ind+2
+ nind=nind+2
+ endif
+ write (iout,*) "i",i," ind",ind," nind",nind
+ enddo
+ if (inct.gt.innt) then
+! DM(ind)=ip4+mp(molnum(inct))/4
+ mnum=molnum(inct)
+ DM(ind)=mp(mnum)/4+ip(mnum)/4
+ if (iabs(itype(inct,mnum)).eq.10.or.molnum(inct).gt.2) then
+ DM(ind)=DM(ind)+msc(itype(inct,molnum(inct)),molnum(inct))
+ ind=ind+1
+ nind=nind+1
+ else
+ mnum=molnum(inct)
+ DM(ind)=DM(ind)+isc(iabs(itype(inct,mnum)),mnum)
+ DM(ind+1)=msc(iabs(itype(inct,mnum)),mnum)+isc(iabs(itype(inct,mnum)),mnum)
+ ind=ind+2
+ nind=nind+2
+ endif
endif
-
-
- ind=1
- do i=nnt,nct
- mnum=molnum(i)
-! if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1 &
-! .and.mnum.eq.5) then
- if (iabs(itype(i,1)).ne.10 .and. &
- iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
- DU1(ind)=-isc(iabs(itype(i,1)),1)
- DU1(ind+1)=0.0d0
- ind=ind+2
- else
- if (mnum.eq.5) then
- DU1(ind)=msc(itype(i,mnum),5))
- else
- DU1(ind)=mp(1)/4-ip4
- ind=ind+1
- endif
- enddo
+ write (iout,*) "ind",ind," nind",nind
+ dimen_chain(ichain)=nind
+ enddo
- ind=1
- do i=nnt,nct-1
- mnum=molnum(i)
-! if (iabs(itype(i,1)).eq.ntyp1) cycle
- write (iout,*) "i",i," itype",itype(i,1),ntyp1
- if (iabs(itype(i,1)).ne.10 .and. &
- iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then
- DU2(ind)=mp(1)/4-ip4
- DU2(ind+1)=0.0d0
- ind=ind+2
- else
- DU2(ind)=0.0d0
- DU2(ind+1)=0.0d0
- ind=ind+1
- endif
- enddo
-#endif
- DMorig=DM
- DU1orig=DU1
- DU2orig=DU2
- write (iout,*) "nind",nind," dimen",dimen
- if (nind.ne.dimen) then
- write (iout,*) "ERROR, dimensions differ"
-#ifdef MPI
- call MPI_Finalize(ierr)
-#endif
- stop
- endif
- write (iout,*) "The upper part of the five-diagonal inertia matrix"
- do i=1,dimen
- if (i.lt.dimen-1) then
- write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
- else if (i.eq.dimen-1) then
- write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
- else
- write (iout,'(i3,3f10.5)') i,DM(i)
- endif
- enddo
+ do ichain=1,nchain
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct
+ mnum=molnum(i)
+ if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ DU1(ind)=-isc(iabs(itype(i,mnum)),mnum)
+ DU1(ind+1)=0.0d0
+ ind=ind+2
+ else
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
+ DU1(ind)=mp(mnum)/4-ip(mnum)/4
+ ind=ind+1
+ endif
+ enddo
+ enddo
- call FDISYP (dimen, DM, DU1, DU2, MARK)
+ do ichain=1,nchain
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct-1
+ mnum=molnum(i)
+! if (iabs(itype(i)).eq.ntyp1) cycle
+!c write (iout,*) "i",i," itype",itype(i),ntyp1
+ if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ DU2(ind)=mp(mnum)/4-ip(mnum)/4
+ DU2(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU2(ind)=0.0d0
+ DU2(ind+1)=0.0d0
+ ind=ind+1
+ endif
+ enddo
+ enddo
+ DMorig=DM
+ DU1orig=DU1
+ DU2orig=DU2
+ if (gmatout) then
+ write (iout,*)"The upper part of the five-diagonal inertia matrix"
+ endif
+ do ichain=1,nchain
+ if (gmatout) write (iout,'(a,i5)') 'Chain',ichain
+ n=dimen_chain(ichain)
+ innt=iposd_chain(ichain)
+ inct=iposd_chain(ichain)+dimen_chain(ichain)-1
+ if (gmatout) then
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
+ endif
+ enddo
+ endif
+ call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),&
+ DU2(innt:inct-1), MARK)
if (mark.eq.-1) then
- write (iout,*) "ERROR: the inertia matrix is not positive definite"
+ write(iout,*)&
+ "ERROR: the inertia matrix is not positive definite for chain",&
+ ichain
#ifdef MPI
call MPI_Finalize(ierr)
#endif
stop
else if (mark.eq.0) then
- write (iout,*) "ERROR: the inertia matrix is singular"
+ write (iout,*)&
+ "ERROR: the inertia matrix is singular for chain",ichain
#ifdef MPI
- call MPI_Finalize(ierr)
+ call MPI_Finalize(ierr)
#endif
else if (mark.eq.1) then
- write (iout,*) "The transformed inertia matrix"
- do i=1,dimen
- if (i.lt.dimen-1) then
- write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
- else if (i.eq.dimen-1) then
- write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+ if (gmatout) then
+ write (iout,*) "The transformed five-diagonal inertia matrix"
+ write (iout,'(a,i5)') 'Chain',ichain
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
else
- write (iout,'(i3,3f10.5)') i,DM(i)
- endif
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
+ endif
enddo
+ endif
endif
+ enddo
! Diagonalization of the pentadiagonal matrix
- allocate(DDU1(2*nres))
- allocate(DDU2(2*nres))
- allocate(DDM(2*nres))
- do i=1,dimen-1
- DDU1(i+1)=DU1orig(i)
- enddo
- do i=1,dimen-2
- DDU2(i+2)=DU2orig(i)
- enddo
- DDM=DMorig
- call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
- if (lprn) then
- write (iout,*) &
- "Eigenvalues of the five-diagonal inertia matrix"
- do i=1,dimen
- write (iout,'(i5,f10.5)') i,geigen(i)
- enddo
- endif
- if (allocated(DDU1)) deallocate(DDU1)
- if (allocated(DDU2)) deallocate(DDU2)
- if (allocated(DDM)) deallocate(DDM)
+#ifdef TIMING
+ time00=MPI_Wtime()
+#endif
+
+
+#endif
+!CRYSTBOND
#else
+ dimen=(nct-nnt+1)+nside
+ dimen1=(nct-nnt)+(nct-nnt+1)
+ dimen3=dimen*3
+ write (iout,*) "nnt",nnt," nct",nct," nside",nside
! Old Gmatrix
#ifdef MPI
if (nfgtasks.gt.1) then
603 FORMAT (I5,4(3F9.3,2x))
604 FORMAT (1H1)
end subroutine MATOUT2
+#ifdef FIVEDIAG
+ subroutine fivediagmult(n,DM,DU1,DU2,x,y)
+ integer n
+ double precision DM(n),DU1(n),DU2(n),x(n),y(n)
+ integer i
+ y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
+ y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
+ do i=3,n-2
+ y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i) &
+ +DU1(i)*x(i+1)+DU2(i)*x(i+2)
+ enddo
+ y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1) &
+ +DU1(n-1)*x(n)
+ y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
+ return
+ end subroutine
+!c---------------------------------------------------------------------------
+ subroutine fivediaginv_mult(ndim,forces,d_a_vec)
+ use energy_data, only:nchain,chain_border,nct,nnt,molnum,&
+ chain_border1,itype
+ integer ndim
+ double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim), &
+ xsolv(ndim),d_a_vec(6*nres)
+ integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev,mnum
+ accel=0.0d0
+ do j=1,3
+!Compute accelerations in Calpha and SC
+ do ichain=1,nchain
+ n=dimen_chain(ichain)
+ iposc=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=iposc,iposc+n-1
+ rs(i-iposc+1)=forces(3*(i-1)+j)
+ enddo
+#ifdef DEBUG
+ write (iout,*) "j",j," chain",ichain
+ write (iout,*) "rs"
+ write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
+ call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
+#ifdef DEBUG
+ write (iout,*) "xsolv"
+ write (iout,'(f10.5)') (xsolv(i),i=1,n)
+#endif
+ ind=1
+ do i=innt,inct
+ mnum=molnum(i)
+ if (itype(i,1).eq.10.or.mnum.gt.2)then
+ accel(j,i)=xsolv(ind)
+ ind=ind+1
+ else
+ accel(j,i)=xsolv(ind)
+ accel(j,i+nres)=xsolv(ind+1)
+ ind=ind+2
+ end if
+ enddo
+ enddo
+ enddo
+!C Convert d_a to virtual-bon-vector basis
+#ifdef DEBUG
+ write (iout,*) "accel in CA-SC basis"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
+ & (accel(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "nnt",nnt
+#endif
+ if (nnt.eq.1) then
+ accel(:,0)=accel(:,1)
+ endif
+ do i=1,nres
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or.mnum.ge.3) then
+ do j=1,3
+ accel(j,i)=accel(j,i+1)-accel(j,i)
+ enddo
+ else
+ do j=1,3
+ accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
+ accel(j,i)=accel(j,i+1)-accel(j,i)
+ enddo
+ end if
+ enddo
+ accel(:,nres)=0.0d0
+ accel(:,nct)=0.0d0
+ accel(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ accel(:,0)=accel(:,1)
+ accel(:,1)=0.0d0
+ endif
+ do ichain=2,nchain
+ accel(:,chain_border1(1,ichain)-1)= &
+ accel(:,chain_border1(1,ichain))
+ accel(:,chain_border1(1,ichain))=0.0d0
+ enddo
+ do ichain=2,nchain
+ accel(:,chain_border1(1,ichain)-1)= &
+ accel(:,chain_border1(1,ichain)-1) &
+ +accel(:,chain_border(2,ichain-1))
+ accel(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+#ifdef DEBUG
+ write (iout,*) "accel in fivediaginv_mult: 1"
+ do i=0,2*nres
+ write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ d_a_vec(j)=accel(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_a_vec(ind+j)=accel(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
+ do j=1,3
+ d_a_vec(ind+j)=accel(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "d_a_vec"
+ write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
+#endif
+ return
+ end subroutine
+
+
+#else
!-----------------------------------------------------------------------------
subroutine ginv_mult(z,d_a_tmp)
return
end subroutine fricmat_mult
!-----------------------------------------------------------------------------
+#endif
!-----------------------------------------------------------------------------
end module REMD
icsa_in=40
!rc for ifc error 118
icsa_pdb=42
+ irotam_end=43
#endif
iscpp=25
icbase=16
real(kind=8),DIMENSION(:),allocatable :: D_ban !(MAXRES6) maxres6=6*maxres
!-----------------------------------------------------------------------------
logical preminim, forceminim ! pre-minimizaation flag
+#ifdef FIVEDIAG
+ integer,parameter :: maxchain=50
+ integer,dimension(maxchain) :: dimen_chain,iposd_chain
+ real(kind=8),dimension(:),allocatable :: DMfric,DU1fric,&
+ DU2fric
+#endif
end module MD_data
integer :: iatsc_s_nucl,iatsc_e_nucl,iatel_s_nucl,iatel_e_nucl,&
iatel_s_vdw_nucl,iatel_e_vdw_nucl,iatscp_s_nucl,iatscp_e_nucl,&
ispp_nucl,iscp_nucl
-
+
! 12/1/95 Array EPS included in the COMMON block.
! common /body/
real(kind=8),dimension(:,:),allocatable :: sigma !(0:ntyp1,0:ntyp1)
integer,dimension(:),allocatable :: nbondterm !(ntyp)
integer,dimension(:),allocatable :: nbondterm_nucl !(ntyp)
+
+
+ integer,dimension(:,:),allocatable :: nterm_scend !(ntyp)
+ real(kind=8),dimension(:,:,:),allocatable:: arotam_end
!-----------------------------------------------------------------------------
! common.local
! Parameters of ab initio-derived potential of virtual-bond-angle bending
real(kind=8),dimension(:,:,:),allocatable :: cref !(3,maxres2+2,maxperm),
real(kind=8),dimension(:,:),allocatable :: crefjlee !(3,maxres2+2),
real(kind=8),dimension(:,:,:),allocatable :: chain_rep !(3,maxres2+2,maxsym)
- integer :: nsup,nstart_sup,nstart_seq,chain_length,iprzes,nperm
+ integer :: nsup,nstart_sup,nstart_seq,iprzes,nperm
integer :: nend_sup,ishift_pdb !wham
real(kind=8) :: rmssing,anatemp !wham
real(kind=8) :: buftubebot, buftubetop,bordtubebot,bordtubetop, &
ibond_nucl,ithep_nucl,irotam_nucl,itorp_nucl,itordp_nucl, &
isidep_nucl,iscpp_nucl,isidep_scbase,isidep_pepbase,isidep_scpho,&
isidep_peppho,iion,iionnucl,iiontran,ilipbond,ilipnonbond,&
- iwaterwater,iwatersc
+ iwaterwater,iwatersc,irotam_end
#ifdef WHAM_RUN
! el wham iounits
integer :: isidep1,ihist,iweight,izsc,idistr
tordname_nucl,sidename_nucl,scpname_nucl, &
sidename_scbase,pepname_pepbase,pepname_scpho,pepname_peppho, &
ionname,ionnuclname,iontranname,lipbondname,lipnonbondname,&
- iwaterwatername,iwaterscname
+ iwaterwatername,iwaterscname,rotname_end
!-----------------------------------------------------------------------
! INP - main input file
! IOUT - list file
! include 'COMMON.VECTORS'
! include 'COMMON.FFIELD'
real(kind=8) :: auxvec(2),auxmat(2,2)
- integer :: i,iti1,iti,k,l
+ integer :: i,iti1,iti,k,l,ii,innt,inct
real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2,cost1,sint1,&
sint1sq,sint1cub,sint1cost1,b1k,b2k,aux
! print *,"in set matrices"
#else
do i=3,nres+1
#endif
+#ifdef FIVEDIAG
+ ii=ireschain(i-2)
+!c write (iout,*) "i",i,i-2," ii",ii
+ if (ii.eq.0) cycle
+ innt=chain_border(1,ii)
+ inct=chain_border(2,ii)
+!c write (iout,*) "i",i,i-2," ii",ii," innt",innt," inct",inct
+!c if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (i.gt. innt+2 .and. i.lt.inct+2) then
+ if (itype(i-2,1).eq.0) then
+ iti = nloctyp
+ else
+ iti = itype2loc(itype(i-2,1))
+ endif
+ else
+ iti=nloctyp
+ endif
+!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. innt+1 .and. i.lt.inct+1) then
+! iti1 = itype2loc(itype(i-1))
+ if (itype(i-1,1).eq.0) then
+ iti1 = nloctyp
+ else
+ iti1 = itype2loc(itype(i-1,1))
+ endif
+ else
+ iti1=nloctyp
+ endif
+#else
if (i.gt. nnt+2 .and. i.lt.nct+2) then
if (itype(i-2,1).eq.0) then
iti = nloctyp
else
iti1=nloctyp
endif
+#endif
! print *,i,itype(i-2,1),iti
#ifdef NEWCORR
cost1=dcos(theta(i-1))
write (iout,*) 'theta=', theta(i-1)
#endif
#else
- if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (i.gt. innt+2 .and. i.lt.inct+2) then
! write(iout,*) "i,",molnum(i),nloctyp
! print *, "i,",molnum(i),i,itype(i-2,1)
if (molnum(i).eq.1) then
#endif
! print *,i,"i"
- if (i .lt. nres+1) then
+ if (i .lt. nres+1 .and. (itype(i-1,1).lt.ntyp1).and.(itype(i-1,1).ne.0)) then
+! if (i .lt. nres+1) then
sin1=dsin(phi(i))
cos1=dcos(phi(i))
sintab(i-2)=sin1
Ug2(2,1,i-2)=0.0d0
Ug2(2,2,i-2)=0.0d0
endif
- if (i .gt. 3 .and. i .lt. nres+1) then
+ if (i .gt. 3) then ! .and. i .lt. nres+1) then
obrot_der(1,i-2)=-sin1
obrot_der(2,i-2)= cos1
Ugder(1,1,i-2)= sin1
do i=ibondp_start,ibondp_end
#ifdef FIVEDIAG
- if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle
+ if (itype(i-1,1).eq.ntyp1 .or. itype(i,1).eq.ntyp1) cycle
diff = vbld(i)-vbldp0
#else
if (itype(i-1,1).eq.ntyp1 .and. itype(i,1).eq.ntyp1) cycle
real(kind=8),dimension(65) :: x
real(kind=8) :: sumene,dsc_i,dp2_i,xx,yy,zz,sumene1,sumene2,sumene3,&
sumene4,s1,s1_6,s2,s2_6,de_dxx,de_dyy,de_dzz,de_dt
- real(kind=8) :: s1_t,s1_6_t,s2_t,s2_6_t
+ real(kind=8) :: s1_t,s1_6_t,s2_t,s2_6_t,gradene
real(kind=8),dimension(3) :: dXX_Ci1,dYY_Ci1,dZZ_Ci1,dXX_Ci,dYY_Ci,&
dZZ_Ci,dXX_XYZ,dYY_XYZ,dZZ_XYZ,dt_dCi,dt_dCi1
!el local variables
- integer :: i,j,k !el,it,nlobit
+ integer :: i,j,k,iti !el,it,nlobit
real(kind=8) :: cosfac2,sinfac2,cosfac,sinfac,escloc,delta
!el real(kind=8) :: time11,time12,time112,theti
!el common /sccalc/ time11,time12,time112,theti,it,nlobit
delta=0.02d0*pi
escloc=0.0D0
do i=loc_start,loc_end
+ gscloc(:,i)=0.0d0
+ gsclocx(:,i)=0.0d0
+! th_gsclocm1(:,i-1)=0.0d0
if (itype(i,1).eq.ntyp1) cycle
costtab(i+1) =dcos(theta(i+1))
sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
it=iabs(itype(i,1))
+ iti=it
+ if (iti.eq.ntyp1 .or. iti.eq.10) cycle
+!c AL 3/30/2022 handle the cases of an isolated-residue chain
+ if (i.eq.nnt .and. itype(i+1,1).eq.ntyp1) cycle
+ if (i.eq.nct .and. itype(i-1,1).eq.ntyp1) cycle
+! costtab(i+1) =dcos(theta(i+1))
if (it.eq.10) goto 1
+#ifdef SC_END
+ if (i.eq.nct .or. itype(i+1,1).eq.ntyp1) then
+!c AL 3/30/2022 handle a sidechain of a loose C-end
+ cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ sumene=arotam_end(0,1,iti)+&
+ tschebyshev(1,nterm_scend(1,iti),arotam_end(1,1,iti),cossc1)
+ escloc=escloc+sumene
+ gradene=gradtschebyshev(0,nterm_scend(1,iti)-1,&
+ arotam_end(1,1,iti),cossc1)
+ gscloc(:,i-1)=gscloc(:,i-1)+&
+ vbld_inv(i)*(dC_norm(:,i+nres)-dC_norm(:,i-1)&
+ *cossc1)*gradene
+ gsclocx(:,i)=gsclocx(:,i)+vbld_inv(i+nres)*&
+ (dC_norm(:,i-1)-dC_norm(:,i+nres)*cossc1)*gradene
+#ifdef ENERGY_DEC
+ if (energy_dec) write (2,'(2hC ,a3,i6,2(a,f10.5))')&
+ restyp(iti,1),i," angle",rad2deg*dacos(cossc1)," escloc",sumene
+#endif
+ else if (i.eq.nnt .or. itype(i-1,1).eq.ntyp1) then
+!c AL 3/30/2022 handle a sidechain of a loose N-end
+ cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ sumene=arotam_end(0,2,iti)+&
+ tschebyshev(1,nterm_scend(2,iti),arotam_end(1,2,iti),cossc)
+ escloc=escloc+sumene
+ gradene=gradtschebyshev(0,nterm_scend(2,iti)-1,&
+ arotam_end(1,2,iti),cossc)
+ gscloc(:,i)=gscloc(:,i)+&
+ vbld_inv(i+1)*(dC_norm(:,i+nres)-dC_norm(:,i)&
+ *cossc)*gradene
+ gsclocx(:,i)=gsclocx(:,i)+vbld_inv(i+nres)*&
+ (dC_norm(:,i)-dC_norm(:,i+nres)*cossc)*gradene
+#ifdef ENERGY_DEC
+ if (energy_dec) write (2,'(2hN ,a3,i6,2(a,f10.5))')
+ & restyp(iti),i," angle",rad2deg*dacos(cossc)," escloc",sumene
+#endif
+ else
+#endif
!
! Compute the axes of tghe local cartesian coordinates system; store in
! x_prime, y_prime and z_prime
! & (gscloc(k,i),k=1,3),(gsclocx(k,i),k=1,3)
! to check gradient call subroutine check_grad
-
+#ifdef SC_END
+ endif
+#endif
1 continue
enddo
return
!d print *,'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
#ifdef FIVEDIAG
call build_fromto(i+1,j+1,fromto)
-c write(iout,'(7hfromto 9f10.5)')((fromto(k,l),l=1,3),k=1,3)
+!c write(iout,'(7hfromto 9f10.5)')((fromto(k,l),l=1,3),k=1,3)
do k=1,3
do l=1,3
tempkl=0.0D0
! call checkintcartgrad
call zerograd
aincr=1.0D-5
- write(iout,*) 'Calling CHECK_ECARTINT.'
+ write(iout,*) 'Calling CHECK_ECARTINT.,kupa'
nf=0
icall=0
call geom_to_var(nvar,x)
call etotal(energia)
etot=energia(0)
call cartgrad
+#ifdef FIVEDIAG
+ call grad_transform
+#endif
icall =1
do i=1,nres
write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
grad_s(j,i)=gcart(j,i)
grad_s(j+3,i)=gxcart(j,i)
write(iout,*) "before movement analytical gradient"
+
+ enddo
+ enddo
do i=1,nres
write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
(gxcart(j,i),j=1,3)
enddo
- enddo
- enddo
else
!- split gradient check
call zerograd
call etotal_long(energia)
!el call enerprint(energia)
call cartgrad
+#ifdef FIVEDIAG
+ call grad_transform
+#endif
icall =1
do i=1,nres
write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
call etotal_short(energia)
call enerprint(energia)
call cartgrad
+#ifdef FIVEDIAG
+ call grad_transform
+#endif
+
icall =1
do i=1,nres
write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
enddo
endif
write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
-! do i=1,nres
+#ifdef FIVEDIAG
+ do i=1,nres
+#else
do i=nnt,nct
+#endif
do j=1,3
if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
call zerograd
call etotal(energia1)
etot1=energia1(0)
- write (iout,*) "ij",i,j," etot1",etot1
+! write (iout,*) "ij",i,j," etot1",etot1
else
!- split gradient
call etotal_long(energia1)
call zerograd
call etotal(energia1)
etot2=energia1(0)
- write (iout,*) "ij",i,j," etot2",etot2
+! write (iout,*) "ij",i,j," etot2",etot2
ggg(j)=(etot1-etot2)/(2*aincr)
else
!- split gradient
#ifdef DEBUG
write (iout,*) "CARGRAD"
#endif
- do i=nres,0,-1
- do j=1,3
- gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+! do i=nres,0,-1
+! do j=1,3
+! gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
- enddo
+! enddo
! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
- enddo
+! enddo
! Correction: dummy residues
- if (nnt.gt.1) then
- do j=1,3
- ! gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
- gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
- enddo
- endif
- if (nct.lt.nres) then
- do j=1,3
- ! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
- gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
- enddo
- endif
+! if (nnt.gt.1) then
+! do j=1,3
+! ! gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
+! gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+! enddo
+! endif
+! if (nct.lt.nres) then
+! do j=1,3
+! ! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+! gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+! enddo
+! endif
+! call grad_transform
#endif
#ifdef TIMING
time_cartgrad=time_cartgrad+MPI_Wtime()-time00
#ifdef MPI
include 'mpif.h'
#endif
- integer i,j,kk
+ integer i,j,kk,mnum
#ifdef DEBUG
write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
write (iout,*) "dC/dX gradient"
enddo
! Correction: dummy residues
do i=2,nres
- if (itype(i-1).eq.ntyp1 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(i)
+ if (itype(i-1,mnum).eq.ntyp1_molec(mnum) .and.&
+ itype(i,mnum).ne.ntyp1_molec(mnum)) then
gcart(:,i)=gcart(:,i)+gcart(:,i-1)
- else if (itype(i-1).ne.ntyp1 .and. itype(i).eq.ntyp1) then
+ else if (itype(i-1,mnum).ne.ntyp1_molec(mnum).and.&
+ itype(i,mnum).eq.ntyp1_molec(mnum)) then
gcart(:,i-1)=gcart(:,i-1)+gcart(:,i)
endif
enddo
-c if (nnt.gt.1) then
-c do j=1,3
-c gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
-c enddo
-c endif
-c if (nct.lt.nres) then
-c do j=1,3
-c! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
-c gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
-c enddo
-c endif
+! if (nnt.gt.1) then
+! do j=1,3
+! gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+! enddo
+! endif
+! if (nct.lt.nres) then
+! do j=1,3
+!! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+! gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+! enddo
+! endif
#ifdef DEBUG
write (iout,*) "CA/SC gradient"
do i=1,nres
#else
do i=3,nres
#endif
- if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1).and.molnum(i).ge.4) then
+ if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1).and.molnum(i).lt.4) then
cost1=dcos(omicron(1,i))
sint1=sqrt(1-cost1*cost1)
cost2=dcos(omicron(2,i))
ctgt1=0.0d0
endif
cosg_inv=1.0d0/cosg
- if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
+! if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
-(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
+(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
- endif
+! endif
! write(iout,*) "just after,close to pi",dphi(j,3,i),&
! sing*(ctgt1*dtheta(j,2,i-1)),ctgt*dtheta(j,1,i), &
! (fac0*vp2(j)+sing*dc_norm(j,i-2)),vbld_inv(i-1)
! Obtaining the gamma derivatives from cosine derivative
else
do j=1,3
- if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
+! if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
dc_norm(j,i-3))/vbld(i-2)
write(iout,*) "just after",dphi(j,3,i),sing,dcosphi(j,3,i)
#endif
!#undef DEBUG
- endif
+! endif
enddo
endif
enddo
use MPI_data
use compare, only:seq_comp,contact
use control
+ implicit none
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
#ifdef MPI
endif
#endif
nct=nres
+ allocate(ireschain(nres))
+ ireschain=0
+ write(iout,*),"before seq2chains",ireschain
+ call seq2chains
+ write(iout,*) "after seq2chains",nchain
+ allocate ( chain_border1(2,nchain))
+ chain_border1(1,1)=1
+ chain_border1(2,1)=chain_border(2,1)+1
+ do i=2,nchain-1
+ chain_border1(1,i)=chain_border(1,i)-1
+ chain_border1(2,i)=chain_border(2,i)+1
+ enddo
+ if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
+ chain_border1(2,nchain)=nres
+ write(iout,*) "nres",nres," nchain",nchain
+ do i=1,nchain
+ write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
+ chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
+ enddo
+ allocate(tabpermchain(50,5040))
+ call chain_symmetry(npermchain,tabpermchain)
print *,'NNT=',NNT,' NCT=',NCT
if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
return
10 stop "Error in fragment file"
end subroutine read_klapaucjusz
+
+!-----------------------------------------------------------
+ subroutine seq2chains
+!c
+!c Split the total UNRES sequence, which has dummy residues separating
+!c the chains, into separate chains. The length of chain ichain is
+!c contained in chain_length(ichain), the first and last non-dummy
+!c residues are in chain_border(1,ichain) and chain_border(2,ichain),
+!c respectively. The lengths pertain to non-dummy residues only.
+!c
+! implicit none
+! include 'DIMENSIONS'
+ use energy_data, only:molnum,nchain,chain_length,ireschain
+ implicit none
+! integer ireschain(nres)
+ integer ii,ichain,i,j,mnum
+ logical new_chain
+ print *,"in seq2"
+ ichain=1
+ new_chain=.true.
+ if (.not.allocated(chain_length)) allocate(chain_length(50))
+ if (.not.allocated(chain_border)) allocate(chain_border(2,50))
+
+ chain_length(ichain)=0
+ ii=1
+ do while (ii.lt.nres)
+ write(iout,*) "in seq2chains",ii,new_chain
+ mnum=molnum(ii)
+ if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
+ if (.not.new_chain) then
+ new_chain=.true.
+ chain_border(2,ichain)=ii-1
+ ichain=ichain+1
+ chain_border(1,ichain)=ii+1
+ chain_length(ichain)=0
+ endif
+ else
+ if (new_chain) then
+ chain_border(1,ichain)=ii
+ new_chain=.false.
+ endif
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ ii=ii+1
+ enddo
+ if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
+ ii=ii-1
+ else
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ if (chain_length(ichain).gt.0) then
+ chain_border(2,ichain)=ii
+ nchain=ichain
+ else
+ nchain=ichain-1
+ endif
+ ireschain=0
+ do i=1,nchain
+ do j=chain_border(1,i),chain_border(2,i)
+ ireschain(j)=i
+ enddo
+ enddo
+ return
+ end subroutine
+!---------------------------------------------------------------------
+ subroutine chain_symmetry(npermchain,tabpermchain)
+!c
+!c Determine chain symmetry. nperm is the number of permutations and
+!c tabperchain contains the allowed permutations of the chains.
+!c
+! implicit none
+! include "DIMENSIONS"
+! include "COMMON.IOUNITS"
+ implicit none
+ integer itemp(50),&
+ npermchain,tabpermchain(50,5040),&
+ tabperm(50,5040),mapchain(50),&
+ iequiv(50,nres),iflag(nres)
+ integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
+ nperm,npermc,ind,mnum
+ if (nchain.eq.1) then
+ npermchain=1
+ tabpermchain(1,1)=1
+!c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
+ return
+ endif
+!c
+!c Look for equivalent chains
+#ifdef DEBUG
+ write(iout,*) "nchain",nchain
+ do i=1,nchain
+ write(iout,*) "chain",i," from",chain_border(1,i),&
+ " to",chain_border(2,i)
+ write(iout,*)&
+ "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
+ enddo
+#endif
+ do i=1,nchain
+ iflag(i)=0
+ enddo
+ nchain_group=0
+ do i=1,nchain
+ if (iflag(i).gt.0) cycle
+ iflag(i)=1
+ nchain_group=nchain_group+1
+ iieq=1
+ iequiv(iieq,nchain_group)=i
+ do j=i+1,nchain
+ if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
+!c k=0
+!c do while(k.lt.chain_length(i) .and.
+!c & itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
+ do k=0,chain_length(i)-1
+!c k=k+1
+ mnum=molnum(k)
+ if (itype(chain_border(1,i)+k,mnum).ne.&
+ itype(chain_border(1,j)+k,mnum)) exit
+ enddo
+ if (k.lt.chain_length(i)) cycle
+ iflag(j)=1
+ iieq=iieq+1
+ iequiv(iieq,nchain_group)=j
+ enddo
+ nequiv(nchain_group)=iieq
+ enddo
+ write(iout,*) "Number of equivalent chain groups:",nchain_group
+ write(iout,*) "Equivalent chain groups"
+ do i=1,nchain_group
+ write(iout,*) "group",i," #members",nequiv(i)," chains",&
+ (iequiv(j,i),j=1,nequiv(i))
+ enddo
+ ind=0
+ do i=1,nchain_group
+ do j=1,nequiv(i)
+ ind=ind+1
+ mapchain(ind)=iequiv(j,i)
+ enddo
+ enddo
+ write (iout,*) "mapchain"
+ do i=1,nchain
+ write (iout,*) i,mapchain(i)
+ enddo
+ ii=0
+ do i=1,nchain_group
+ call permut(nequiv(i),nperm,tabperm)
+ if (ii.eq.0) then
+ ii=nequiv(i)
+ npermchain=nperm
+ do j=1,nperm
+ do k=1,ii
+ tabpermchain(k,j)=iequiv(tabperm(k,j),i)
+ enddo
+ enddo
+ else
+ npermc=npermchain
+ npermchain=npermchain*nperm
+ ind=0
+ do k=1,nperm
+ do j=1,npermc
+ ind=ind+1
+ do l=1,ii
+ tabpermchain(l,ind)=tabpermchain(l,j)
+ enddo
+ do l=1,nequiv(i)
+ tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
+ enddo
+ enddo
+ enddo
+ ii=ii+nequiv(i)
+ endif
+ enddo
+ do i=1,npermchain
+ do j=1,nchain
+ itemp(mapchain(j))=tabpermchain(j,i)
+ enddo
+ do j=1,nchain
+ tabpermchain(j,i)=itemp(j)
+ enddo
+ enddo
+ write(iout,*) "Number of chain permutations",npermchain
+ write(iout,*) "Permutations"
+ do i=1,npermchain
+ write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
+ enddo
+ return
+ end
+!c---------------------------------------------------------------------
+ integer function tperm(i,iperm,tabpermchain)
+! implicit none
+! include 'DIMENSIONS'
+ integer i,iperm
+ integer tabpermchain(50,5040)
+ if (i.eq.0) then
+ tperm=0
+ else
+ tperm=tabpermchain(i,iperm)
+ endif
+ return
+ end function
+
!-----------------------------------------------------------------------------
end module io
!-----------------------------------------------------------------------------
! permut.F
!-----------------------------------------------------------------------------
- subroutine permut(isym)
-
- use geometry_data, only: tabperm
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
-! include 'COMMON.LOCAL'
-! include 'COMMON.VAR'
-! include 'COMMON.CHAIN'
-! include 'COMMON.INTERACT'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.GEO'
-! include 'COMMON.NAMES'
-! include 'COMMON.CONTROL'
-
- integer :: n,isym
+ subroutine permut(isym,nperm,tabperm)
+!c integer maxperm,maxsym
+!c parameter (maxperm=3628800)
+!c parameter (maxsym=10)
+! include "DIMENSIONS"
+ integer n,a,tabperm,nperm,kkk,i,isym
! logical nextp
-!el external nextp
- integer,dimension(isym) :: a
-! parameter(n=symetr)
-!el local variables
- integer :: kkk,i
-
+! external nextp
+ dimension a(isym),tabperm(50,5040)
n=isym
+ nperm=1
if (n.eq.1) then
tabperm(1,1)=1
return
endif
+ do i=2,n
+ nperm=nperm*i
+ enddo
kkk=0
do i=1,n
a(i)=i
enddo
- 10 print *,(a(i),i=1,n)
+ 10 continue
+!c print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n)
kkk=kkk+1
do i=1,n
- tabperm(kkk,i)=a(i)
-! write (iout,*) "tututu", kkk
+ tabperm(i,kkk)=a(i)
enddo
if(nextp(n,a)) go to 10
return
- end subroutine permut
+ end subroutine
+
+
!-----------------------------------------------------------------------------
logical function nextp(n,a)
dsc_inv(i)=1.0D0/dsc(i)
endif
enddo
+! ip(1)=0.0001d0
+! isc(:,1)=0.0001d0
#endif
read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
do i=1,ntyp_molec(2)
enddo
endif
enddo
+#ifdef SC_END
+ allocate(nterm_scend(2,ntyp))
+ allocate(arotam_end(0:6,2,ntyp))
+ nterm_scend=0
+ arotam_end=0.0d0
+ read (irotam_end,*) ijunk
+!c write (iout,*) "ijunk",ijunk
+ do i=1,ntyp
+ if (i.eq.10) cycle
+ do j=1,2
+ read (irotam_end,'(a)')
+ read (irotam_end,*) nterm_scend(j,i)
+!c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
+ do k=0,nterm_scend(j,i)
+ read (irotam_end,*) ijunk,arotam_end(k,j,i)
+!c write (iout,*) "k",k," arotam",arotam_end(k,j,i)
+ enddo
+ enddo
+ enddo
+!c lprint=.true.
+ if (lprint) then
+ write (iout,'(a)') &
+ "Parameters of the local potentials of sidechain ends"
+ do i=1,ntyp
+ write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
+ restyp(i,1),restyp(i,1)
+ do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
+ write (iout,'(i5,2f20.10)') &
+ j,arotam_end(j,1,i),arotam_end(j,2,i)
+ enddo
+ enddo
+ endif
+!c lprint=.false.
+#endif
+
!---------reading nucleic acid parameters for rotamers-------------------
allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
do i=1,ntyp_molec(2)
sigiso2(i,j)=sigiso2(j,i)
! print *,"ATU",sigma(j,i),sigma(i,j),i,j
nstate(i,j) = nstate(j,i)
- dtail(1,i,j) = dtail(1,j,i)
- dtail(2,i,j) = dtail(2,j,i)
+ dtail(1,i,j) = dtail(2,j,i)
+ dtail(2,i,j) = dtail(1,j,i)
DO k = 1, 4
alphasur(k,i,j) = alphasur(k,j,i)
wstate(k,i,j) = wstate(k,j,i)
cou=1
write (iout,*) "symetr", symetr
do i=1,nres
- lll=lll+1
+ lll=lll+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
- if (i.gt.1) then
- if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
- chain_length=lll-1
- kkk=kkk+1
+! if (i.gt.1) then
+! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
+! chain_length=lll-1
+! kkk=kkk+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
- lll=1
- endif
- endif
+! lll=1
+! endif
+! endif
do j=1,3
cref(j,i,cou)=c(j,i)
cref(j,i+nres,cou)=c(j,i+nres)
endif
enddo
enddo
- write (iout,*) chain_length
- if (chain_length.eq.0) chain_length=nres
- do j=1,3
- chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
- chain_rep(j,chain_length+nres,symetr) &
- =chain_rep(j,chain_length+nres,1)
- enddo
! diagnostic
! write (iout,*) "spraw lancuchy",chain_length,symetr
! do i=1,4
dc(j,0)=c(j,1)
enddo
- if (symetr.gt.1) then
- call permut(symetr)
- nperm=1
- do i=1,symetr
- nperm=nperm*i
- enddo
- do i=1,nperm
- write(iout,*) (tabperm(i,kkk),kkk=1,4)
- enddo
- do i=1,nperm
- cou=0
- do kkk=1,symetr
- icha=tabperm(i,kkk)
- write (iout,*) i,icha
- do lll=1,chain_length
- cou=cou+1
- if (cou.le.nres) then
- do j=1,3
- kupa=mod(lll,chain_length)
- iprzes=(kkk-1)*chain_length+lll
- if (kupa.eq.0) kupa=chain_length
- write (iout,*) "kupa", kupa
- cref(j,iprzes,i)=chain_rep(j,kupa,icha)
- cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
- enddo
- endif
- enddo
- enddo
- enddo
- endif
+! if (symetr.gt.1) then
+! call permut(symetr)
+! nperm=1
+! do i=1,symetr
+! nperm=nperm*i
+! enddo
+! do i=1,nperm
+! write(iout,*) (tabperm(i,kkk),kkk=1,4)
+! enddo
+! do i=1,nperm
+! cou=0
+! do kkk=1,symetr
+! icha=tabperm(i,kkk)
+! write (iout,*) i,icha
+! do lll=1,chain_length
+! cou=cou+1
+! if (cou.le.nres) then
+! do j=1,3
+! kupa=mod(lll,chain_length)
+! iprzes=(kkk-1)*chain_length+lll
+! if (kupa.eq.0) kupa=chain_length
+! write (iout,*) "kupa", kupa
+! cref(j,iprzes,i)=chain_rep(j,kupa,icha)
+! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
+! enddo
+! endif
+! enddo
+! enddo
+! enddo
+! endif
!-koniec robienia kopii
! diag
do kkk=1,nperm
! print *,"Processor",myrank," opened file ITHEP"
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
! print *,"Processor",myrank," opened file IROTAM"
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
open (ithep,file=thetname,status='old')
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
open (ithep,file=thetname,status='old',action='read')
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
call getenv_loc('TORDPAR',tordname)
use MD_data
use energy
use MDyn, only:setup_fricmat
+#ifndef FIVEDIAG
use REMD, only:fricmat_mult,ginv_mult
+#endif
#ifdef MPI
include "mpif.h"
#endif
! write (2,*) "After sum_gradient"
! write (2,*) "dimen",dimen," dimen3",dimen3
! call flush(2)
+#ifndef FIVEDIAG
else if (iorder.eq.4) then
call ginv_mult(z,d_a_tmp)
else if (iorder.eq.5) then
! call flush(2)
! write (iout,*) "My chunk of ginv_block"
! call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
+#endif
else if (iorder.eq.6) then
call int_from_cart1(.false.)
else if (iorder.eq.7) then
call chainbuild_cart
else if (iorder.eq.8) then
call intcartderiv
+#ifndef FIVEDIAG
else if (iorder.eq.9) then
call fricmat_mult(z,d_a_tmp)
+#endif
else if (iorder.eq.10) then
call setup_fricmat
endif
projects / unres4.git / commitdiff