8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm(ntyp)) !(ntyp)
780 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
782 allocate(msc(ntyp+1)) !(ntyp+1)
783 allocate(isc(ntyp+1)) !(ntyp+1)
784 allocate(restok(ntyp+1)) !(ntyp+1)
785 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(long_r_sidechain(ntyp))
787 allocate(short_r_sidechain(ntyp))
792 read (ibond,*) vbldp0,akp,mp,ip,pstok
795 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
796 dsc(i) = vbldsc0(1,i)
800 dsc_inv(i)=1.0D0/dsc(i)
804 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
806 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
807 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
808 dsc(i) = vbldsc0(1,i)
812 dsc_inv(i)=1.0D0/dsc(i)
817 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
818 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
820 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
822 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
823 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
825 write (iout,'(13x,3f10.5)') &
826 vbldsc0(j,i),aksc(j,i),abond0(j,i)
830 !----------------------------------------------------
831 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
832 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
833 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
834 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
835 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
836 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
846 allocate(liptranene(ntyp))
847 !C reading lipid parameters
848 write (iout,*) "iliptranpar",iliptranpar
850 read(iliptranpar,*) pepliptran
853 read(iliptranpar,*) liptranene(i)
854 print *,liptranene(i)
860 ! Read the parameters of the probability distribution/energy expression
861 ! of the virtual-bond valence angles theta
864 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
865 (bthet(j,i,1,1),j=1,2)
866 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
867 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
868 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
872 athet(1,i,1,-1)=athet(1,i,1,1)
873 athet(2,i,1,-1)=athet(2,i,1,1)
874 bthet(1,i,1,-1)=-bthet(1,i,1,1)
875 bthet(2,i,1,-1)=-bthet(2,i,1,1)
876 athet(1,i,-1,1)=-athet(1,i,1,1)
877 athet(2,i,-1,1)=-athet(2,i,1,1)
878 bthet(1,i,-1,1)=bthet(1,i,1,1)
879 bthet(2,i,-1,1)=bthet(2,i,1,1)
883 athet(1,i,-1,-1)=athet(1,-i,1,1)
884 athet(2,i,-1,-1)=-athet(2,-i,1,1)
885 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
886 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
887 athet(1,i,-1,1)=athet(1,-i,1,1)
888 athet(2,i,-1,1)=-athet(2,-i,1,1)
889 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
890 bthet(2,i,-1,1)=bthet(2,-i,1,1)
891 athet(1,i,1,-1)=-athet(1,-i,1,1)
892 athet(2,i,1,-1)=athet(2,-i,1,1)
893 bthet(1,i,1,-1)=bthet(1,-i,1,1)
894 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
899 polthet(j,i)=polthet(j,-i)
902 gthet(j,i)=gthet(j,-i)
910 'Parameters of the virtual-bond valence angles:'
911 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
912 ' ATHETA0 ',' A1 ',' A2 ',&
915 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
916 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
918 write (iout,'(/a/9x,5a/79(1h-))') &
919 'Parameters of the expression for sigma(theta_c):',&
920 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
921 ' ALPH3 ',' SIGMA0C '
923 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
924 (polthet(j,i),j=0,3),sigc0(i)
926 write (iout,'(/a/9x,5a/79(1h-))') &
927 'Parameters of the second gaussian:',&
928 ' THETA0 ',' SIGMA0 ',' G1 ',&
931 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
932 sig0(i),(gthet(j,i),j=1,3)
936 'Parameters of the virtual-bond valence angles:'
937 write (iout,'(/a/9x,5a/79(1h-))') &
938 'Coefficients of expansion',&
939 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
940 ' b1*10^1 ',' b2*10^1 '
942 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
943 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
944 (10*bthet(j,i,1,1),j=1,2)
946 write (iout,'(/a/9x,5a/79(1h-))') &
947 'Parameters of the expression for sigma(theta_c):',&
948 ' alpha0 ',' alph1 ',' alph2 ',&
949 ' alhp3 ',' sigma0c '
951 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
952 (polthet(j,i),j=0,3),sigc0(i)
954 write (iout,'(/a/9x,5a/79(1h-))') &
955 'Parameters of the second gaussian:',&
956 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
959 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
960 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
966 ! Read the parameters of Utheta determined from ab initio surfaces
967 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
969 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
970 ntheterm3,nsingle,ndouble
971 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
973 !----------------------------------------------------
974 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
975 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
976 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
977 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
978 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
979 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
980 !(maxtheterm,-maxthetyp1:maxthetyp1,&
981 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
982 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
983 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
984 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
985 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
986 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
987 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
988 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
989 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
990 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
991 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
992 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
993 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
994 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
995 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
996 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
997 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
999 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1001 ithetyp(i)=-ithetyp(-i)
1004 aa0thet(:,:,:,:)=0.0d0
1005 aathet(:,:,:,:,:)=0.0d0
1006 bbthet(:,:,:,:,:,:)=0.0d0
1007 ccthet(:,:,:,:,:,:)=0.0d0
1008 ddthet(:,:,:,:,:,:)=0.0d0
1009 eethet(:,:,:,:,:,:)=0.0d0
1010 ffthet(:,:,:,:,:,:,:)=0.0d0
1011 ggthet(:,:,:,:,:,:,:)=0.0d0
1013 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1015 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1016 ! VAR:1=non-glicyne non-proline 2=proline
1017 ! VAR:negative values for D-aminoacid
1019 do j=-nthetyp,nthetyp
1020 do k=-nthetyp,nthetyp
1021 read (ithep,'(6a)',end=111,err=111) res1
1022 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1023 ! VAR: aa0thet is variable describing the average value of Foureir
1024 ! VAR: expansion series
1025 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1026 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1027 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1028 read (ithep,*,end=111,err=111) &
1029 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1030 read (ithep,*,end=111,err=111) &
1031 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1032 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1033 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1034 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1036 read (ithep,*,end=111,err=111) &
1037 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1038 ffthet(lll,llll,ll,i,j,k,iblock),&
1039 ggthet(llll,lll,ll,i,j,k,iblock),&
1040 ggthet(lll,llll,ll,i,j,k,iblock),&
1041 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1046 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1047 ! coefficients of theta-and-gamma-dependent terms are zero.
1048 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1049 ! RECOMENTDED AFTER VERSION 3.3)
1053 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1054 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1056 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1057 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1060 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1062 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1065 ! AND COMMENT THE LOOPS BELOW
1069 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1070 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1072 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1073 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1076 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1078 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1083 ! Substitution for D aminoacids from symmetry.
1086 do j=-nthetyp,nthetyp
1087 do k=-nthetyp,nthetyp
1088 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1090 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1094 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1095 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1096 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1097 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1103 ffthet(llll,lll,ll,i,j,k,iblock)= &
1104 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1105 ffthet(lll,llll,ll,i,j,k,iblock)= &
1106 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1107 ggthet(llll,lll,ll,i,j,k,iblock)= &
1108 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1109 ggthet(lll,llll,ll,i,j,k,iblock)= &
1110 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1119 ! Control printout of the coefficients of virtual-bond-angle potentials
1122 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1127 write (iout,'(//4a)') &
1128 'Type ',onelett(i),onelett(j),onelett(k)
1129 write (iout,'(//a,10x,a)') " l","a[l]"
1130 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1131 write (iout,'(i2,1pe15.5)') &
1132 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1134 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1135 "b",l,"c",l,"d",l,"e",l
1137 write (iout,'(i2,4(1pe15.5))') m,&
1138 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1139 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1143 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1144 "f+",l,"f-",l,"g+",l,"g-",l
1147 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1148 ffthet(n,m,l,i,j,k,iblock),&
1149 ffthet(m,n,l,i,j,k,iblock),&
1150 ggthet(n,m,l,i,j,k,iblock),&
1151 ggthet(m,n,l,i,j,k,iblock)
1161 write (2,*) "Start reading THETA_PDB",ithep_pdb
1163 ! write (2,*) 'i=',i
1164 read (ithep_pdb,*,err=111,end=111) &
1165 a0thet(i),(athet(j,i,1,1),j=1,2),&
1166 (bthet(j,i,1,1),j=1,2)
1167 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1168 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1169 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1170 sigc0(i)=sigc0(i)**2
1173 athet(1,i,1,-1)=athet(1,i,1,1)
1174 athet(2,i,1,-1)=athet(2,i,1,1)
1175 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1176 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1177 athet(1,i,-1,1)=-athet(1,i,1,1)
1178 athet(2,i,-1,1)=-athet(2,i,1,1)
1179 bthet(1,i,-1,1)=bthet(1,i,1,1)
1180 bthet(2,i,-1,1)=bthet(2,i,1,1)
1183 a0thet(i)=a0thet(-i)
1184 athet(1,i,-1,-1)=athet(1,-i,1,1)
1185 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1186 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1187 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1188 athet(1,i,-1,1)=athet(1,-i,1,1)
1189 athet(2,i,-1,1)=-athet(2,-i,1,1)
1190 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1191 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1192 athet(1,i,1,-1)=-athet(1,-i,1,1)
1193 athet(2,i,1,-1)=athet(2,-i,1,1)
1194 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1195 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1196 theta0(i)=theta0(-i)
1200 polthet(j,i)=polthet(j,-i)
1203 gthet(j,i)=gthet(j,-i)
1206 write (2,*) "End reading THETA_PDB"
1211 !-------------------------------------------
1212 allocate(nlob(ntyp1)) !(ntyp1)
1213 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1214 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1215 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1222 gaussc(:,:,:,:)=0.0D0
1226 ! Read the parameters of the probability distribution/energy expression
1227 ! of the side chains.
1230 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1234 dsc_inv(i)=1.0D0/dsc(i)
1245 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1246 ((blower(k,l,1),l=1,k),k=1,3)
1247 censc(1,1,-i)=censc(1,1,i)
1248 censc(2,1,-i)=censc(2,1,i)
1249 censc(3,1,-i)=-censc(3,1,i)
1251 read (irotam,*,end=112,err=112) bsc(j,i)
1252 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1253 ((blower(k,l,j),l=1,k),k=1,3)
1254 censc(1,j,-i)=censc(1,j,i)
1255 censc(2,j,-i)=censc(2,j,i)
1256 censc(3,j,-i)=-censc(3,j,i)
1257 ! BSC is amplitude of Gaussian
1264 akl=akl+blower(k,m,j)*blower(l,m,j)
1268 if (((k.eq.3).and.(l.ne.3)) &
1269 .or.((l.eq.3).and.(k.ne.3))) then
1270 gaussc(k,l,j,-i)=-akl
1271 gaussc(l,k,j,-i)=-akl
1273 gaussc(k,l,j,-i)=akl
1274 gaussc(l,k,j,-i)=akl
1283 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1286 if (nlobi.gt.0) then
1288 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1289 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1290 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1291 'log h',(bsc(j,i),j=1,nlobi)
1292 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1293 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1295 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1296 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1299 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1300 write (iout,'(a,f10.4,4(16x,f10.4))') &
1301 'Center ',(bsc(j,i),j=1,nlobi)
1302 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1311 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1312 ! added by Urszula Kozlowska 07/11/2007
1314 !el Maximum number of SC local term fitting function coefficiants
1315 !el integer,parameter :: maxsccoef=65
1317 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1320 read (irotam,*,end=112,err=112)
1322 read (irotam,*,end=112,err=112)
1325 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1330 ! Read the parameters of the probability distribution/energy expression
1331 ! of the side chains.
1333 write (2,*) "Start reading ROTAM_PDB"
1335 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1339 dsc_inv(i)=1.0D0/dsc(i)
1350 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1351 ((blower(k,l,1),l=1,k),k=1,3)
1353 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1354 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1355 ((blower(k,l,j),l=1,k),k=1,3)
1362 akl=akl+blower(k,m,j)*blower(l,m,j)
1372 write (2,*) "End reading ROTAM_PDB"
1378 ! Read torsional parameters in old format
1380 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1382 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1383 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1384 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1386 !el from energy module--------
1387 allocate(v1(nterm_old,ntortyp,ntortyp))
1388 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1389 !el---------------------------
1394 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1400 write (iout,'(/a/)') 'Torsional constants:'
1403 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1404 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1410 ! Read torsional parameters
1412 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1413 read (itorp,*,end=113,err=113) ntortyp
1414 !el from energy module---------
1415 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1416 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1418 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1419 allocate(vlor2(maxlor,ntortyp,ntortyp))
1420 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1421 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1423 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1424 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1425 !el---------------------------
1428 !el---------------------------
1430 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1432 itortyp(i)=-itortyp(-i)
1434 itortyp(ntyp1)=ntortyp
1435 itortyp(-ntyp1)=-ntortyp
1437 write (iout,*) 'ntortyp',ntortyp
1439 do j=-ntortyp+1,ntortyp-1
1440 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1442 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1443 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1446 do k=1,nterm(i,j,iblock)
1447 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1449 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1450 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1451 v0ij=v0ij+si*v1(k,i,j,iblock)
1453 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1454 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1455 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1457 do k=1,nlor(i,j,iblock)
1458 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1459 vlor2(k,i,j),vlor3(k,i,j)
1460 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1463 v0(-i,-j,iblock)=v0ij
1469 write (iout,'(/a/)') 'Torsional constants:'
1471 do i=-ntortyp,ntortyp
1472 do j=-ntortyp,ntortyp
1473 write (iout,*) 'ityp',i,' jtyp',j
1474 write (iout,*) 'Fourier constants'
1475 do k=1,nterm(i,j,iblock)
1476 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1479 write (iout,*) 'Lorenz constants'
1480 do k=1,nlor(i,j,iblock)
1481 write (iout,'(3(1pe15.5))') &
1482 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1488 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1490 ! 6/23/01 Read parameters for double torsionals
1492 !el from energy module------------
1493 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1494 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1495 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1496 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1497 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1498 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1499 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1500 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1501 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1502 !---------------------------------
1506 do j=-ntortyp+1,ntortyp-1
1507 do k=-ntortyp+1,ntortyp-1
1508 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1509 ! write (iout,*) "OK onelett",
1512 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1513 .or. t3.ne.toronelet(k)) then
1514 write (iout,*) "Error in double torsional parameter file",&
1517 call MPI_Finalize(Ierror)
1519 stop "Error in double torsional parameter file"
1521 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1522 ntermd_2(i,j,k,iblock)
1523 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1524 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1525 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1526 ntermd_1(i,j,k,iblock))
1527 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1528 ntermd_1(i,j,k,iblock))
1529 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1530 ntermd_1(i,j,k,iblock))
1531 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1532 ntermd_1(i,j,k,iblock))
1533 ! Martix of D parameters for one dimesional foureir series
1534 do l=1,ntermd_1(i,j,k,iblock)
1535 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1536 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1537 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1538 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1539 ! write(iout,*) "whcodze" ,
1540 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1542 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1543 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1544 v2s(m,l,i,j,k,iblock),&
1545 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1546 ! Martix of D parameters for two dimesional fourier series
1547 do l=1,ntermd_2(i,j,k,iblock)
1549 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1550 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1551 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1552 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1561 write (iout,*) 'Constants for double torsionals'
1564 do j=-ntortyp+1,ntortyp-1
1565 do k=-ntortyp+1,ntortyp-1
1566 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1567 ' nsingle',ntermd_1(i,j,k,iblock),&
1568 ' ndouble',ntermd_2(i,j,k,iblock)
1570 write (iout,*) 'Single angles:'
1571 do l=1,ntermd_1(i,j,k,iblock)
1572 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1573 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1574 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1575 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1578 write (iout,*) 'Pairs of angles:'
1579 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1580 do l=1,ntermd_2(i,j,k,iblock)
1581 write (iout,'(i5,20f10.5)') &
1582 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1585 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1586 do l=1,ntermd_2(i,j,k,iblock)
1587 write (iout,'(i5,20f10.5)') &
1588 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1589 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1598 ! Read of Side-chain backbone correlation parameters
1599 ! Modified 11 May 2012 by Adasko
1602 read (isccor,*,end=119,err=119) nsccortyp
1604 !el from module energy-------------
1605 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1606 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1607 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1608 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1609 !-----------------------------------
1611 !el from module energy-------------
1612 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1614 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1616 isccortyp(i)=-isccortyp(-i)
1618 iscprol=isccortyp(20)
1619 ! write (iout,*) 'ntortyp',ntortyp
1621 !c maxinter is maximum interaction sites
1622 !el from module energy---------
1623 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1624 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1625 -nsccortyp:nsccortyp))
1626 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1627 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1628 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1629 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1630 !-----------------------------------
1634 read (isccor,*,end=119,err=119) &
1635 nterm_sccor(i,j),nlor_sccor(i,j)
1641 nterm_sccor(-i,j)=nterm_sccor(i,j)
1642 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1643 nterm_sccor(i,-j)=nterm_sccor(i,j)
1644 do k=1,nterm_sccor(i,j)
1645 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1647 if (j.eq.iscprol) then
1648 if (i.eq.isccortyp(10)) then
1649 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1650 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1652 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1653 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1654 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1655 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1656 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1657 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1658 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1659 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1662 if (i.eq.isccortyp(10)) then
1663 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1664 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1666 if (j.eq.isccortyp(10)) then
1667 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1668 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1670 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1671 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1672 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1673 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1674 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1675 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1679 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1680 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1681 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1682 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1685 do k=1,nlor_sccor(i,j)
1686 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1687 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1688 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1689 (1+vlor3sccor(k,i,j)**2)
1691 v0sccor(l,i,j)=v0ijsccor
1692 v0sccor(l,-i,j)=v0ijsccor1
1693 v0sccor(l,i,-j)=v0ijsccor2
1694 v0sccor(l,-i,-j)=v0ijsccor3
1700 !el from module energy-------------
1701 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1703 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1704 ! write (iout,*) 'ntortyp',ntortyp
1706 !c maxinter is maximum interaction sites
1707 !el from module energy---------
1708 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1709 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1710 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1711 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1712 !-----------------------------------
1716 read (isccor,*,end=119,err=119) &
1717 nterm_sccor(i,j),nlor_sccor(i,j)
1721 do k=1,nterm_sccor(i,j)
1722 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1724 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1727 do k=1,nlor_sccor(i,j)
1728 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1729 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1730 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1731 (1+vlor3sccor(k,i,j)**2)
1733 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1741 write (iout,'(/a/)') 'Torsional constants:'
1744 write (iout,*) 'ityp',i,' jtyp',j
1745 write (iout,*) 'Fourier constants'
1746 do k=1,nterm_sccor(i,j)
1747 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1749 write (iout,*) 'Lorenz constants'
1750 do k=1,nlor_sccor(i,j)
1751 write (iout,'(3(1pe15.5))') &
1752 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1759 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1760 ! interaction energy of the Gly, Ala, and Pro prototypes.
1765 write (iout,*) "Coefficients of the cumulants"
1767 read (ifourier,*) nloctyp
1768 !write(iout,*) "nloctyp",nloctyp
1769 !el from module energy-------
1770 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1771 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1772 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1773 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1774 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1775 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1776 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1777 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1781 !--------------------------------
1784 read (ifourier,*,end=115,err=115)
1785 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1787 write (iout,*) 'Type',i
1788 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1798 B1tilde(1,-i) =-b(3)
1800 ! b1tilde(1,i)=0.0d0
1801 ! b1tilde(2,i)=0.0d0
1826 Ctilde(1,2,-i)=-b(9)
1830 ! Ctilde(1,1,i)=0.0d0
1831 ! Ctilde(1,2,i)=0.0d0
1832 ! Ctilde(2,1,i)=0.0d0
1833 ! Ctilde(2,2,i)=0.0d0
1851 Dtilde(1,2,-i)=-b(8)
1855 ! Dtilde(1,1,i)=0.0d0
1856 ! Dtilde(1,2,i)=0.0d0
1857 ! Dtilde(2,1,i)=0.0d0
1858 ! Dtilde(2,2,i)=0.0d0
1859 EE(1,1,i)= b(10)+b(11)
1860 EE(2,2,i)=-b(10)+b(11)
1861 EE(2,1,i)= b(12)-b(13)
1862 EE(1,2,i)= b(12)+b(13)
1863 EE(1,1,-i)= b(10)+b(11)
1864 EE(2,2,-i)=-b(10)+b(11)
1865 EE(2,1,-i)=-b(12)+b(13)
1866 EE(1,2,-i)=-b(12)-b(13)
1872 ! ee(2,1,i)=ee(1,2,i)
1876 write (iout,*) 'Type',i
1878 write(iout,*) B1(1,i),B1(2,i)
1880 write(iout,*) B2(1,i),B2(2,i)
1883 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1887 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1891 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1896 ! Read electrostatic-interaction parameters
1901 write (iout,'(/a)') 'Electrostatic interaction constants:'
1902 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1903 'IT','JT','APP','BPP','AEL6','AEL3'
1905 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1906 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1907 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1908 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1913 app (i,j)=epp(i,j)*rri*rri
1914 bpp (i,j)=-2.0D0*epp(i,j)*rri
1915 ael6(i,j)=elpp6(i,j)*4.2D0**6
1916 ael3(i,j)=elpp3(i,j)*4.2D0**3
1918 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1924 ! Read side-chain interaction parameters.
1926 !el from module energy - COMMON.INTERACT-------
1927 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1928 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1929 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1930 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1931 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1932 allocate(epslip(ntyp,ntyp))
1940 !--------------------------------
1942 read (isidep,*,end=117,err=117) ipot,expon
1943 if (ipot.lt.1 .or. ipot.gt.5) then
1944 write (iout,'(2a)') 'Error while reading SC interaction',&
1945 'potential file - unknown potential type.'
1947 call MPI_Finalize(Ierror)
1953 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1954 ', exponents are ',expon,2*expon
1955 ! goto (10,20,30,30,40) ipot
1957 !----------------------- LJ potential ---------------------------------
1959 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1960 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961 (sigma0(i),i=1,ntyp)
1963 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1964 write (iout,'(a/)') 'The epsilon array:'
1965 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1966 write (iout,'(/a)') 'One-body parameters:'
1967 write (iout,'(a,4x,a)') 'residue','sigma'
1968 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1971 !----------------------- LJK potential --------------------------------
1973 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1974 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1975 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1977 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1978 write (iout,'(a/)') 'The epsilon array:'
1979 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1980 write (iout,'(/a)') 'One-body parameters:'
1981 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1982 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1986 !---------------------- GB or BP potential -----------------------------
1990 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1992 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1993 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1994 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1995 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1997 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2000 ! For the GB potential convert sigma'**2 into chi'
2003 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2007 write (iout,'(/a/)') 'Parameters of the BP potential:'
2008 write (iout,'(a/)') 'The epsilon array:'
2009 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010 write (iout,'(/a)') 'One-body parameters:'
2011 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2013 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
2014 chip(i),alp(i),i=1,ntyp)
2017 !--------------------- GBV potential -----------------------------------
2019 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2020 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2021 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2022 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2024 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2025 write (iout,'(a/)') 'The epsilon array:'
2026 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2027 write (iout,'(/a)') 'One-body parameters:'
2028 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2029 's||/s_|_^2',' chip ',' alph '
2030 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2031 sigii(i),chip(i),alp(i),i=1,ntyp)
2034 write(iout,*)"Wrong ipot"
2039 !-----------------------------------------------------------------------
2040 ! Calculate the "working" parameters of SC interactions.
2042 !el from module energy - COMMON.INTERACT-------
2043 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2044 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2045 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2046 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2048 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2061 sc_aa_tube_par(:)=0.0d0
2062 sc_bb_tube_par(:)=0.0d0
2064 !--------------------------------
2069 epslip(i,j)=epslip(j,i)
2074 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2075 sigma(j,i)=sigma(i,j)
2076 rs0(i,j)=dwa16*sigma(i,j)
2080 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2081 'Working parameters of the SC interactions:',&
2082 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2087 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2096 sigeps=dsign(1.0D0,epsij)
2098 aa_aq(i,j)=epsij*rrij*rrij
2099 bb_aq(i,j)=-sigeps*epsij*rrij
2100 aa_aq(j,i)=aa_aq(i,j)
2101 bb_aq(j,i)=bb_aq(i,j)
2102 epsijlip=epslip(i,j)
2103 sigeps=dsign(1.0D0,epsijlip)
2104 epsijlip=dabs(epsijlip)
2105 aa_lip(i,j)=epsijlip*rrij*rrij
2106 bb_lip(i,j)=-sigeps*epsijlip*rrij
2107 aa_lip(j,i)=aa_lip(i,j)
2108 bb_lip(j,i)=bb_lip(i,j)
2109 !C write(iout,*) aa_lip
2111 sigt1sq=sigma0(i)**2
2112 sigt2sq=sigma0(j)**2
2115 ratsig1=sigt2sq/sigt1sq
2116 ratsig2=1.0D0/ratsig1
2117 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2118 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2119 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2123 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2124 sigmaii(i,j)=rsum_max
2125 sigmaii(j,i)=rsum_max
2127 ! sigmaii(i,j)=r0(i,j)
2128 ! sigmaii(j,i)=r0(i,j)
2130 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2131 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2132 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2133 augm(i,j)=epsij*r_augm**(2*expon)
2134 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2141 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2142 restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2143 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2147 write(iout,*) "tube param"
2148 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2149 ccavtubpep,dcavtubpep,tubetranenepep
2150 sigmapeptube=sigmapeptube**6
2151 sigeps=dsign(1.0D0,epspeptube)
2152 epspeptube=dabs(epspeptube)
2153 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2154 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2155 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2157 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2158 ccavtub(i),dcavtub(i),tubetranene(i)
2159 sigmasctube=sigmasctube**6
2160 sigeps=dsign(1.0D0,epssctube)
2161 epssctube=dabs(epssctube)
2162 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2163 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2164 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2167 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2172 ! Define the SC-p interaction constants (hard-coded; old style)
2175 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2177 ! aad(i,1)=0.3D0*4.0D0**12
2178 ! Following line for constants currently implemented
2179 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2180 aad(i,1)=1.5D0*4.0D0**12
2181 ! aad(i,1)=0.17D0*5.6D0**12
2183 ! "Soft" SC-p repulsion
2185 ! Following line for constants currently implemented
2186 ! aad(i,1)=0.3D0*4.0D0**6
2187 ! "Hard" SC-p repulsion
2188 bad(i,1)=3.0D0*4.0D0**6
2189 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2198 ! 8/9/01 Read the SC-p interaction constants from file
2201 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2204 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2205 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2206 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2207 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2211 write (iout,*) "Parameters of SC-p interactions:"
2213 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2214 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2220 ! Define the constants of the disulfide bridge
2224 ! Old arbitrary potential - commented out.
2229 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2230 ! energy surface of diethyl disulfide.
2231 ! A. Liwo and U. Kozlowska, 11/24/03
2248 write (iout,'(/a)') "Disulfide bridge parameters:"
2249 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2250 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2251 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2252 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2256 111 write (iout,*) "Error reading bending energy parameters."
2258 112 write (iout,*) "Error reading rotamer energy parameters."
2260 113 write (iout,*) "Error reading torsional energy parameters."
2262 114 write (iout,*) "Error reading double torsional energy parameters."
2264 115 write (iout,*) &
2265 "Error reading cumulant (multibody energy) parameters."
2267 116 write (iout,*) "Error reading electrostatic energy parameters."
2269 117 write (iout,*) "Error reading side chain interaction parameters."
2271 118 write (iout,*) "Error reading SCp interaction parameters."
2273 119 write (iout,*) "Error reading SCCOR parameters"
2276 call MPI_Finalize(Ierror)
2280 end subroutine parmread
2282 !-----------------------------------------------------------------------------
2284 !-----------------------------------------------------------------------------
2285 subroutine printmat(ldim,m,n,iout,key,a)
2288 character(len=3),dimension(n) :: key
2289 real(kind=8),dimension(ldim,n) :: a
2291 integer :: i,j,k,m,iout,nlim
2295 write (iout,1000) (key(k),k=i,nlim)
2297 1000 format (/5x,8(6x,a3))
2298 1020 format (/80(1h-)/)
2300 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2303 1010 format (a3,2x,8(f9.4))
2305 end subroutine printmat
2306 !-----------------------------------------------------------------------------
2308 !-----------------------------------------------------------------------------
2310 ! Read the PDB file and convert the peptide geometry into virtual-chain
2313 use energy_data, only: itype
2317 use control, only: rescode
2318 ! implicit real*8 (a-h,o-z)
2319 ! include 'DIMENSIONS'
2320 ! include 'COMMON.LOCAL'
2321 ! include 'COMMON.VAR'
2322 ! include 'COMMON.CHAIN'
2323 ! include 'COMMON.INTERACT'
2324 ! include 'COMMON.IOUNITS'
2325 ! include 'COMMON.GEO'
2326 ! include 'COMMON.NAMES'
2327 ! include 'COMMON.CONTROL'
2328 ! include 'COMMON.DISTFIT'
2329 ! include 'COMMON.SETUP'
2330 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2332 logical :: lprn=.true.,fail
2333 real(kind=8),dimension(3) :: e1,e2,e3
2334 real(kind=8) :: dcj,efree_temp
2335 character(len=3) :: seq,res
2336 character(len=5) :: atom
2337 character(len=80) :: card
2338 real(kind=8),dimension(3,20) :: sccor
2339 integer :: kkk,lll,icha,kupa !rescode,
2342 integer,dimension(2,maxres/3) :: hfrag_alloc
2343 integer,dimension(4,maxres/3) :: bfrag_alloc
2344 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2350 ! write (2,*) "UNRES_PDB",unres_pdb
2358 !-----------------------------
2359 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2360 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2363 read (ipdbin,'(a80)',end=10) card
2364 ! write (iout,'(a)') card
2365 if (card(:5).eq.'HELIX') then
2368 read(card(22:25),*) hfrag(1,nhfrag)
2369 read(card(34:37),*) hfrag(2,nhfrag)
2371 if (card(:5).eq.'SHEET') then
2374 read(card(24:26),*) bfrag(1,nbfrag)
2375 read(card(35:37),*) bfrag(2,nbfrag)
2376 !rc----------------------------------------
2377 !rc to be corrected !!!
2378 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2379 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2380 !rc----------------------------------------
2382 if (card(:3).eq.'END') then
2384 else if (card(:3).eq.'TER') then
2388 itype(ires_old)=ntyp1
2389 itype(ires_old-1)=ntyp1
2391 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2394 dc(j,ires)=sccor(j,iii)
2397 call sccenter(ires,iii,sccor)
2403 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2404 ! Fish out the ATOM cards.
2405 if (index(card(1:4),'ATOM').gt.0) then
2406 read (card(12:16),*) atom
2407 ! write (iout,*) "! ",atom," !",ires
2408 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2409 read (card(23:26),*) ires
2410 read (card(18:20),'(a3)') res
2411 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2412 ! & " ires_old",ires_old
2413 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2414 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2415 if (ires-ishift+ishift1.ne.ires_old) then
2416 ! Calculate the CM of the preceding residue.
2417 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2419 ! write (iout,*) "Calculating sidechain center iii",iii
2422 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2425 call sccenter(ires_old,iii,sccor)
2429 ! Start new residue.
2430 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2433 else if (ibeg.eq.1) then
2434 write (iout,*) "BEG ires",ires
2436 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2440 ires=ires-ishift+ishift1
2442 ! write (iout,*) "ishift",ishift," ires",ires,&
2443 ! " ires_old",ires_old
2445 else if (ibeg.eq.2) then
2447 ishift=-ires_old+ires-1 !!!!!
2448 ishift1=ishift1-1 !!!!!
2449 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2450 ires=ires-ishift+ishift1
2454 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2455 ires=ires-ishift+ishift1
2458 if (res.eq.'ACE' .or. res.eq.'NHE') then
2461 itype(ires)=rescode(ires,res,0)
2464 ires=ires-ishift+ishift1
2466 ! write (iout,*) "ires_old",ires_old," ires",ires
2467 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2470 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2471 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2472 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2473 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2474 ! write (iout,*) "backbone ",atom
2476 write (iout,'(2i3,2x,a,3f8.3)') &
2477 ires,itype(ires),res,(c(j,ires),j=1,3)
2481 sccor(j,iii)=c(j,ires)
2483 ! write (*,*) card(23:27),ires,itype(ires)
2484 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2485 atom.ne.'N' .and. atom.ne.'C' .and. &
2486 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2487 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2488 ! write (iout,*) "sidechain ",atom
2490 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2494 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2495 if (ires.eq.0) return
2496 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2500 ! write (iout,*) i,itype(i)
2501 ! if (itype(i).eq.ntyp1) then
2502 ! write (iout,*) "dummy",i,itype(i)
2504 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2505 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2509 if (itype(i).eq.ntyp1) then
2510 if (itype(i+1).eq.ntyp1) then
2511 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2512 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
2513 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
2515 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2516 ! print *,i,'tu dochodze'
2517 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2525 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2529 dcj=(c(j,i-2)-c(j,i-3))/2.0
2530 if (dcj.eq.0) dcj=1.23591524223
2535 else !itype(i+1).eq.ntyp1
2537 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2538 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2545 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2549 dcj=(c(j,i+3)-c(j,i+2))/2.0
2550 if (dcj.eq.0) dcj=1.23591524223
2555 endif !itype(i+1).eq.ntyp1
2556 endif !itype.eq.ntyp1
2559 ! Calculate the CM of the last side chain.
2563 dc(j,ires)=sccor(j,iii)
2566 call sccenter(ires,iii,sccor)
2572 if (itype(nres).ne.10) then
2576 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2577 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2584 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2588 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2589 c(j,nres)=c(j,nres-1)+dcj
2590 c(j,2*nres)=c(j,nres)
2594 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2596 if (nres.ne.nres0) then
2597 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2599 stop "Error nres value in WHAM input"
2602 !---------------------------------
2603 !el reallocate tables
2606 ! hfrag_alloc(j,i)=hfrag(j,i)
2609 ! bfrag_alloc(j,i)=bfrag(j,i)
2615 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2616 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2617 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2618 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2622 ! hfrag(j,i)=hfrag_alloc(j,i)
2627 ! bfrag(j,i)=bfrag_alloc(j,i)
2630 !el end reallocate tables
2631 !---------------------------------
2639 c(j,2*nres)=c(j,nres)
2641 if (itype(1).eq.ntyp1) then
2645 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2646 call refsys(2,3,4,e1,e2,e3,fail)
2653 c(j,1)=c(j,2)-1.9d0*e2(j)
2657 dcj=(c(j,4)-c(j,3))/2.0
2663 ! Copy the coordinates to reference coordinates
2669 ! Calculate internal coordinates.
2671 write (iout,'(/a)') &
2672 "Cartesian coordinates of the reference structure"
2673 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2674 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2676 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2677 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2678 (c(j,ires+nres),j=1,3)
2681 ! znamy już nres wiec mozna alokowac tablice
2682 ! Calculate internal coordinates.
2683 if(me.eq.king.or..not.out1file)then
2684 write (iout,'(a)') &
2685 "Backbone and SC coordinates as read from the PDB"
2687 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2688 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2689 (c(j,nres+ires),j=1,3)
2693 if(.not.allocated(vbld)) then
2694 allocate(vbld(2*nres))
2699 if(.not.allocated(vbld_inv)) then
2700 allocate(vbld_inv(2*nres))
2706 if(.not.allocated(theta)) then
2707 allocate(theta(nres+2))
2711 if(.not.allocated(phi)) allocate(phi(nres+2))
2712 if(.not.allocated(alph)) allocate(alph(nres+2))
2713 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2714 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2715 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2716 if(.not.allocated(costtab)) allocate(costtab(nres))
2717 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2718 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2719 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2720 if(.not.allocated(xxref)) allocate(xxref(nres))
2721 if(.not.allocated(yyref)) allocate(yyref(nres))
2722 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2723 if(.not.allocated(dc_norm)) then
2724 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2725 allocate(dc_norm(3,0:2*nres+2))
2729 call int_from_cart(.true.,.false.)
2730 call sc_loc_geom(.false.)
2732 thetaref(i)=theta(i)
2742 dc(j,i)=c(j,i+1)-c(j,i)
2743 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2748 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2749 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2751 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2755 ! Copy the coordinates to reference coordinates
2756 ! Splits to single chain if occurs
2760 ! cref(j,i,cou)=c(j,i)
2764 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2765 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2766 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2767 !-----------------------------
2773 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2775 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2778 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2783 cref(j,i,cou)=c(j,i)
2784 cref(j,i+nres,cou)=c(j,i+nres)
2786 chain_rep(j,lll,kkk)=c(j,i)
2787 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2791 write (iout,*) chain_length
2792 if (chain_length.eq.0) chain_length=nres
2794 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2795 chain_rep(j,chain_length+nres,symetr) &
2796 =chain_rep(j,chain_length+nres,1)
2799 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2801 ! do kkk=1,chain_length
2802 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2806 ! makes copy of chains
2807 write (iout,*) "symetr", symetr
2812 if (symetr.gt.1) then
2819 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2825 ! write (iout,*) i,icha
2826 do lll=1,chain_length
2828 if (cou.le.nres) then
2830 kupa=mod(lll,chain_length)
2831 iprzes=(kkk-1)*chain_length+lll
2832 if (kupa.eq.0) kupa=chain_length
2833 ! write (iout,*) "kupa", kupa
2834 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2835 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2842 !-koniec robienia kopii
2845 write (iout,*) "nowa struktura", nperm
2847 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2849 cref(3,i,kkk),cref(1,nres+i,kkk),&
2850 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2852 100 format (//' alpha-carbon coordinates ',&
2853 ' centroid coordinates'/ &
2854 ' ', 6X,'X',11X,'Y',11X,'Z', &
2855 10X,'X',11X,'Y',11X,'Z')
2856 110 format (a,'(',i3,')',6f12.5)
2862 bfrag(i,j)=bfrag(i,j)-ishift
2868 hfrag(i,j)=hfrag(i,j)-ishift
2874 end subroutine readpdb
2875 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2876 !-----------------------------------------------------------------------------
2878 !-----------------------------------------------------------------------------
2879 subroutine read_control
2893 use random, only: random_init
2894 ! implicit real*8 (a-h,o-z)
2895 ! include 'DIMENSIONS'
2897 use prng, only:prng_restart
2899 logical :: OKRandom!, prng_restart
2902 ! include 'COMMON.IOUNITS'
2903 ! include 'COMMON.TIME1'
2904 ! include 'COMMON.THREAD'
2905 ! include 'COMMON.SBRIDGE'
2906 ! include 'COMMON.CONTROL'
2907 ! include 'COMMON.MCM'
2908 ! include 'COMMON.MAP'
2909 ! include 'COMMON.HEADER'
2910 ! include 'COMMON.CSA'
2911 ! include 'COMMON.CHAIN'
2912 ! include 'COMMON.MUCA'
2913 ! include 'COMMON.MD'
2914 ! include 'COMMON.FFIELD'
2915 ! include 'COMMON.INTERACT'
2916 ! include 'COMMON.SETUP'
2917 !el integer :: KDIAG,ICORFL,IXDR
2918 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2919 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2920 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2921 ! character(len=80) :: ucase
2922 character(len=640) :: controlcard
2924 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2930 read (INP,'(a)') titel
2931 call card_concat(controlcard,.true.)
2932 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2933 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2934 call reada(controlcard,'SEED',seed,0.0D0)
2935 call random_init(seed)
2936 ! Set up the time limit (caution! The time must be input in minutes!)
2937 read_cart=index(controlcard,'READ_CART').gt.0
2938 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2939 call readi(controlcard,'SYM',symetr,1)
2940 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2941 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2942 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2943 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2944 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2945 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2946 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2947 call reada(controlcard,'DRMS',drms,0.1D0)
2948 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2949 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2950 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2951 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2952 write (iout,'(a,f10.1)')'DRMS = ',drms
2953 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2954 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2956 call readi(controlcard,'NZ_START',nz_start,0)
2957 call readi(controlcard,'NZ_END',nz_end,0)
2958 call readi(controlcard,'IZ_SC',iz_sc,0)
2959 timlim=60.0D0*timlim
2960 safety = 60.0d0*safety
2963 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2964 !C SHIELD keyword sets if the shielding effect of side-chains is used
2965 !C 0 denots no shielding is used all peptide are equally despite the
2966 !C solvent accesible area
2967 !C 1 the newly introduced function
2968 !C 2 reseved for further possible developement
2969 call readi(controlcard,'SHIELD',shield_mode,0)
2970 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2971 write(iout,*) "shield_mode",shield_mode
2972 !C Varibles set size of box
2973 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
2974 write (iout,*) "with_theta_constr ",with_theta_constr
2975 AFMlog=(index(controlcard,'AFM'))
2976 selfguide=(index(controlcard,'SELFGUIDE'))
2977 print *,'AFMlog',AFMlog,selfguide,"KUPA"
2978 call readi(controlcard,'GENCONSTR',genconstr,0)
2979 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2980 call reada(controlcard,'BOXY',boxysize,100.0d0)
2981 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2982 call readi(controlcard,'TUBEMOD',tubemode,0)
2983 write (iout,*) TUBEmode,"TUBEMODE"
2984 if (TUBEmode.gt.0) then
2985 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2986 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2987 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2988 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2989 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2990 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2991 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2992 buftubebot=bordtubebot+tubebufthick
2993 buftubetop=bordtubetop-tubebufthick
2996 ! CUTOFFF ON ELECTROSTATICS
2997 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2998 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2999 write(iout,*) "R_CUT_ELE=",r_cut_ele
3000 ! Lipidic parameters
3001 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3002 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3003 if (lipthick.gt.0.0d0) then
3004 bordliptop=(boxzsize+lipthick)/2.0
3005 bordlipbot=bordliptop-lipthick
3006 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3007 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3008 buflipbot=bordlipbot+lipbufthick
3009 bufliptop=bordliptop-lipbufthick
3010 if ((lipbufthick*2.0d0).gt.lipthick) &
3011 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3012 endif !lipthick.gt.0
3013 write(iout,*) "bordliptop=",bordliptop
3014 write(iout,*) "bordlipbot=",bordlipbot
3015 write(iout,*) "bufliptop=",bufliptop
3016 write(iout,*) "buflipbot=",buflipbot
3017 write (iout,*) "SHIELD MODE",shield_mode
3019 !C-------------------------
3020 minim=(index(controlcard,'MINIMIZE').gt.0)
3021 dccart=(index(controlcard,'CART').gt.0)
3022 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3023 overlapsc=.not.overlapsc
3024 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3025 searchsc=.not.searchsc
3026 sideadd=(index(controlcard,'SIDEADD').gt.0)
3027 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3028 outpdb=(index(controlcard,'PDBOUT').gt.0)
3029 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3030 pdbref=(index(controlcard,'PDBREF').gt.0)
3031 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3032 indpdb=index(controlcard,'PDBSTART')
3033 extconf=(index(controlcard,'EXTCONF').gt.0)
3034 call readi(controlcard,'IPRINT',iprint,0)
3035 call readi(controlcard,'MAXGEN',maxgen,10000)
3036 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3037 call readi(controlcard,"KDIAG",kdiag,0)
3038 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3039 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3040 write (iout,*) "RESCALE_MODE",rescale_mode
3041 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3042 if (index(controlcard,'REGULAR').gt.0.0D0) then
3043 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3047 if (index(controlcard,'CHECKGRAD').gt.0) then
3049 if (index(controlcard,'CART').gt.0) then
3051 elseif (index(controlcard,'CARINT').gt.0) then
3056 elseif (index(controlcard,'THREAD').gt.0) then
3058 call readi(controlcard,'THREAD',nthread,0)
3059 if (nthread.gt.0) then
3060 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3063 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3064 stop 'Error termination in Read_Control.'
3066 else if (index(controlcard,'MCMA').gt.0) then
3068 else if (index(controlcard,'MCEE').gt.0) then
3070 else if (index(controlcard,'MULTCONF').gt.0) then
3072 else if (index(controlcard,'MAP').gt.0) then
3074 call readi(controlcard,'MAP',nmap,0)
3075 else if (index(controlcard,'CSA').gt.0) then
3077 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3079 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3082 !fcm else if (index(controlcard,'MCMF').gt.0) then
3084 else if (index(controlcard,'SOFTREG').gt.0) then
3086 else if (index(controlcard,'CHECK_BOND').gt.0) then
3088 else if (index(controlcard,'TEST').gt.0) then
3090 else if (index(controlcard,'MD').gt.0) then
3092 else if (index(controlcard,'RE ').gt.0) then
3096 lmuca=index(controlcard,'MUCA').gt.0
3097 call readi(controlcard,'MUCADYN',mucadyn,0)
3098 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3099 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3101 write (iout,*) 'MUCADYN=',mucadyn
3102 write (iout,*) 'MUCASMOOTH=',muca_smooth
3105 iscode=index(controlcard,'ONE_LETTER')
3106 indphi=index(controlcard,'PHI')
3107 indback=index(controlcard,'BACK')
3108 iranconf=index(controlcard,'RAND_CONF')
3109 i2ndstr=index(controlcard,'USE_SEC_PRED')
3110 gradout=index(controlcard,'GRADOUT').gt.0
3111 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3112 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3113 if (me.eq.king .or. .not.out1file ) &
3114 write (iout,*) "DISTCHAINMAX",distchainmax
3116 if(me.eq.king.or..not.out1file) &
3117 write (iout,'(2a)') diagmeth(kdiag),&
3118 ' routine used to diagonalize matrices.'
3119 if (shield_mode.gt.0) then
3121 !C VSolvSphere the volume of solving sphere
3123 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3124 !C there will be no distinction between proline peptide group and normal peptide
3125 !C group in case of shielding parameters
3126 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3127 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3128 write (iout,*) VSolvSphere,VSolvSphere_div
3129 !C long axis of side chain
3131 long_r_sidechain(i)=vbldsc0(1,i)
3132 short_r_sidechain(i)=sigma0(i)
3133 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3139 end subroutine read_control
3140 !-----------------------------------------------------------------------------
3141 subroutine read_REMDpar
3143 ! Read REMD settings
3150 use control_data, only:out1file
3151 ! implicit real*8 (a-h,o-z)
3152 ! include 'DIMENSIONS'
3153 ! include 'COMMON.IOUNITS'
3154 ! include 'COMMON.TIME1'
3155 ! include 'COMMON.MD'
3158 !el include 'COMMON.LANGEVIN'
3160 !el include 'COMMON.LANGEVIN.lang0'
3162 ! include 'COMMON.INTERACT'
3163 ! include 'COMMON.NAMES'
3164 ! include 'COMMON.GEO'
3165 ! include 'COMMON.REMD'
3166 ! include 'COMMON.CONTROL'
3167 ! include 'COMMON.SETUP'
3168 ! character(len=80) :: ucase
3169 character(len=320) :: controlcard
3170 character(len=3200) :: controlcard1
3171 integer :: iremd_m_total
3174 ! real(kind=8) :: var,ene
3176 if(me.eq.king.or..not.out1file) &
3177 write (iout,*) "REMD setup"
3179 call card_concat(controlcard,.true.)
3180 call readi(controlcard,"NREP",nrep,3)
3181 call readi(controlcard,"NSTEX",nstex,1000)
3182 call reada(controlcard,"RETMIN",retmin,10.0d0)
3183 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3184 mremdsync=(index(controlcard,'SYNC').gt.0)
3185 call readi(controlcard,"NSYN",i_sync_step,100)
3186 restart1file=(index(controlcard,'REST1FILE').gt.0)
3187 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3188 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3189 if(max_cache_traj_use.gt.max_cache_traj) &
3190 max_cache_traj_use=max_cache_traj
3191 if(me.eq.king.or..not.out1file) then
3192 !d if (traj1file) then
3193 !rc caching is in testing - NTWX is not ignored
3194 !d write (iout,*) "NTWX value is ignored"
3195 !d write (iout,*) " trajectory is stored to one file by master"
3196 !d write (iout,*) " before exchange at NSTEX intervals"
3198 write (iout,*) "NREP= ",nrep
3199 write (iout,*) "NSTEX= ",nstex
3200 write (iout,*) "SYNC= ",mremdsync
3201 write (iout,*) "NSYN= ",i_sync_step
3202 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3205 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3206 if (index(controlcard,'TLIST').gt.0) then
3208 call card_concat(controlcard1,.true.)
3209 read(controlcard1,*) (remd_t(i),i=1,nrep)
3210 if(me.eq.king.or..not.out1file) &
3211 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3214 if (index(controlcard,'MLIST').gt.0) then
3216 call card_concat(controlcard1,.true.)
3217 read(controlcard1,*) (remd_m(i),i=1,nrep)
3218 if(me.eq.king.or..not.out1file) then
3219 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3222 iremd_m_total=iremd_m_total+remd_m(i)
3224 write (iout,*) 'Total number of replicas ',iremd_m_total
3227 if(me.eq.king.or..not.out1file) &
3228 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3230 end subroutine read_REMDpar
3231 !-----------------------------------------------------------------------------
3232 subroutine read_MDpar
3236 use control_data, only: r_cut,rlamb,out1file
3238 use geometry_data, only: pi
3240 ! implicit real*8 (a-h,o-z)
3241 ! include 'DIMENSIONS'
3242 ! include 'COMMON.IOUNITS'
3243 ! include 'COMMON.TIME1'
3244 ! include 'COMMON.MD'
3247 !el include 'COMMON.LANGEVIN'
3249 !el include 'COMMON.LANGEVIN.lang0'
3251 ! include 'COMMON.INTERACT'
3252 ! include 'COMMON.NAMES'
3253 ! include 'COMMON.GEO'
3254 ! include 'COMMON.SETUP'
3255 ! include 'COMMON.CONTROL'
3256 ! include 'COMMON.SPLITELE'
3257 ! character(len=80) :: ucase
3258 character(len=320) :: controlcard
3263 call card_concat(controlcard,.true.)
3264 call readi(controlcard,"NSTEP",n_timestep,1000000)
3265 call readi(controlcard,"NTWE",ntwe,100)
3266 call readi(controlcard,"NTWX",ntwx,1000)
3267 call reada(controlcard,"DT",d_time,1.0d-1)
3268 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3269 call reada(controlcard,"DAMAX",damax,1.0d1)
3270 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3271 call readi(controlcard,"LANG",lang,0)
3272 RESPA = index(controlcard,"RESPA") .gt. 0
3273 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3274 ntime_split0=ntime_split
3275 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3276 ntime_split0=ntime_split
3277 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3278 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3279 rest = index(controlcard,"REST").gt.0
3280 tbf = index(controlcard,"TBF").gt.0
3281 usampl = index(controlcard,"USAMPL").gt.0
3282 mdpdb = index(controlcard,"MDPDB").gt.0
3283 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3284 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3285 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3286 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3287 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3288 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3289 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3290 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3291 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3292 large = index(controlcard,"LARGE").gt.0
3293 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3294 rattle = index(controlcard,"RATTLE").gt.0
3295 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3301 if(me.eq.king.or..not.out1file) then
3303 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3305 write (iout,'(a)') "The units are:"
3306 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3307 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3308 " acceleration: angstrom/(48.9 fs)**2"
3309 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3311 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3312 write (iout,'(a60,f10.5,a)') &
3313 "Initial time step of numerical integration:",d_time,&
3315 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3317 write (iout,'(2a,i4,a)') &
3318 "A-MTS algorithm used; initial time step for fast-varying",&
3319 " short-range forces split into",ntime_split," steps."
3320 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3321 r_cut," lambda",rlamb
3323 write (iout,'(2a,f10.5)') &
3324 "Maximum acceleration threshold to reduce the time step",&
3325 "/increase split number:",damax
3326 write (iout,'(2a,f10.5)') &
3327 "Maximum predicted energy drift to reduce the timestep",&
3328 "/increase split number:",edriftmax
3329 write (iout,'(a60,f10.5)') &
3330 "Maximum velocity threshold to reduce velocities:",dvmax
3331 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3332 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3333 if (rattle) write (iout,'(a60)') &
3334 "Rattle algorithm used to constrain the virtual bonds"
3338 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3339 call reada(controlcard,"RWAT",rwat,1.4d0)
3340 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3341 surfarea=index(controlcard,"SURFAREA").gt.0
3342 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3343 if(me.eq.king.or..not.out1file)then
3344 write (iout,'(/a,$)') "Langevin dynamics calculation"
3346 write (iout,'(a/)') &
3347 " with direct integration of Langevin equations"
3348 else if (lang.eq.2) then
3349 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3350 else if (lang.eq.3) then
3351 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3352 else if (lang.eq.4) then
3353 write (iout,'(a/)') " in overdamped mode"
3355 write (iout,'(//a,i5)') &
3356 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3359 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3360 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3361 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3362 write (iout,'(a60,f10.5)') &
3363 "Scaling factor of the friction forces:",scal_fric
3364 if (surfarea) write (iout,'(2a,i10,a)') &
3365 "Friction coefficients will be scaled by solvent-accessible",&
3366 " surface area every",reset_fricmat," steps."
3368 ! Calculate friction coefficients and bounds of stochastic forces
3369 eta=6*pi*cPoise*etawat
3370 if(me.eq.king.or..not.out1file) &
3371 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3373 gamp=scal_fric*(pstok+rwat)*eta
3374 stdfp=dsqrt(2*Rb*t_bath/d_time)
3375 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3377 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3378 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3380 if(me.eq.king.or..not.out1file)then
3381 write (iout,'(/2a/)') &
3382 "Radii of site types and friction coefficients and std's of",&
3383 " stochastic forces of fully exposed sites"
3384 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3386 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3387 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3391 if(me.eq.king.or..not.out1file)then
3392 write (iout,'(a)') "Berendsen bath calculation"
3393 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3394 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3396 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3397 count_reset_moment," steps"
3399 write (iout,'(a,i10,a)') &
3400 "Velocities will be reset at random every",count_reset_vel,&
3404 if(me.eq.king.or..not.out1file) &
3405 write (iout,'(a31)') "Microcanonical mode calculation"
3407 if(me.eq.king.or..not.out1file)then
3408 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3410 write(iout,*) "MD running with constraints."
3411 write(iout,*) "Equilibration time ", eq_time, " mtus."
3412 write(iout,*) "Constraining ", nfrag," fragments."
3413 write(iout,*) "Length of each fragment, weight and q0:"
3415 write (iout,*) "Set of restraints #",iset
3417 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3418 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3420 write(iout,*) "constraints between ", npair, "fragments."
3421 write(iout,*) "constraint pairs, weights and q0:"
3423 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3424 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3426 write(iout,*) "angle constraints within ", nfrag_back,&
3427 "backbone fragments."
3428 write(iout,*) "fragment, weights:"
3430 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3431 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3432 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3435 iset=mod(kolor,nset)+1
3438 if(me.eq.king.or..not.out1file) &
3439 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3441 end subroutine read_MDpar
3442 !-----------------------------------------------------------------------------
3446 ! implicit real*8 (a-h,o-z)
3447 ! include 'DIMENSIONS'
3448 ! include 'COMMON.MAP'
3449 ! include 'COMMON.IOUNITS'
3450 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3451 character(len=80) :: mapcard !,ucase
3454 ! real(kind=8) :: var,ene
3457 read (inp,'(a)') mapcard
3458 mapcard=ucase(mapcard)
3459 if (index(mapcard,'PHI').gt.0) then
3461 else if (index(mapcard,'THE').gt.0) then
3463 else if (index(mapcard,'ALP').gt.0) then
3465 else if (index(mapcard,'OME').gt.0) then
3468 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3469 stop 'Error - illegal variable spec in MAP card.'
3471 call readi (mapcard,'RES1',res1(imap),0)
3472 call readi (mapcard,'RES2',res2(imap),0)
3473 if (res1(imap).eq.0) then
3474 res1(imap)=res2(imap)
3475 else if (res2(imap).eq.0) then
3476 res2(imap)=res1(imap)
3478 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3479 write (iout,'(a)') &
3480 'Error - illegal definition of variable group in MAP.'
3481 stop 'Error - illegal definition of variable group in MAP.'
3483 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3484 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3485 call readi(mapcard,'NSTEP',nstep(imap),0)
3486 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3487 write (iout,'(a)') &
3488 'Illegal boundary and/or step size specification in MAP.'
3489 stop 'Illegal boundary and/or step size specification in MAP.'
3493 end subroutine map_read
3494 !-----------------------------------------------------------------------------
3497 use control_data, only: vdisulf
3499 ! implicit real*8 (a-h,o-z)
3500 ! include 'DIMENSIONS'
3501 ! include 'COMMON.IOUNITS'
3502 ! include 'COMMON.GEO'
3503 ! include 'COMMON.CSA'
3504 ! include 'COMMON.BANK'
3505 ! include 'COMMON.CONTROL'
3506 ! character(len=80) :: ucase
3507 character(len=620) :: mcmcard
3509 ! integer :: ntf,ik,iw_pdb
3510 ! real(kind=8) :: var,ene
3512 call card_concat(mcmcard,.true.)
3514 call readi(mcmcard,'NCONF',nconf,50)
3515 call readi(mcmcard,'NADD',nadd,0)
3516 call readi(mcmcard,'JSTART',jstart,1)
3517 call readi(mcmcard,'JEND',jend,1)
3518 call readi(mcmcard,'NSTMAX',nstmax,500000)
3519 call readi(mcmcard,'N0',n0,1)
3520 call readi(mcmcard,'N1',n1,6)
3521 call readi(mcmcard,'N2',n2,4)
3522 call readi(mcmcard,'N3',n3,0)
3523 call readi(mcmcard,'N4',n4,0)
3524 call readi(mcmcard,'N5',n5,0)
3525 call readi(mcmcard,'N6',n6,10)
3526 call readi(mcmcard,'N7',n7,0)
3527 call readi(mcmcard,'N8',n8,0)
3528 call readi(mcmcard,'N9',n9,0)
3529 call readi(mcmcard,'N14',n14,0)
3530 call readi(mcmcard,'N15',n15,0)
3531 call readi(mcmcard,'N16',n16,0)
3532 call readi(mcmcard,'N17',n17,0)
3533 call readi(mcmcard,'N18',n18,0)
3535 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3537 call readi(mcmcard,'NDIFF',ndiff,2)
3538 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3539 call readi(mcmcard,'IS1',is1,1)
3540 call readi(mcmcard,'IS2',is2,8)
3541 call readi(mcmcard,'NRAN0',nran0,4)
3542 call readi(mcmcard,'NRAN1',nran1,2)
3543 call readi(mcmcard,'IRR',irr,1)
3544 call readi(mcmcard,'NSEED',nseed,20)
3545 call readi(mcmcard,'NTOTAL',ntotal,10000)
3546 call reada(mcmcard,'CUT1',cut1,2.0d0)
3547 call reada(mcmcard,'CUT2',cut2,5.0d0)
3548 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3549 call readi(mcmcard,'ICMAX',icmax,3)
3550 call readi(mcmcard,'IRESTART',irestart,0)
3551 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3554 call reada(mcmcard,'DELE',dele,20.0d0)
3555 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3556 call readi(mcmcard,'IREF',iref,0)
3557 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3558 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3559 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3560 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3561 write (iout,*) "NCONF_IN",nconf_in
3563 end subroutine csaread
3564 !-----------------------------------------------------------------------------
3568 use control_data, only: MaxMoveType
3571 ! implicit real*8 (a-h,o-z)
3572 ! include 'DIMENSIONS'
3573 ! include 'COMMON.MCM'
3574 ! include 'COMMON.MCE'
3575 ! include 'COMMON.IOUNITS'
3576 ! character(len=80) :: ucase
3577 character(len=320) :: mcmcard
3580 ! real(kind=8) :: var,ene
3582 call card_concat(mcmcard,.true.)
3583 call readi(mcmcard,'MAXACC',maxacc,100)
3584 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3585 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3586 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3587 call readi(mcmcard,'MAXREPM',maxrepm,200)
3588 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3589 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3590 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3591 call reada(mcmcard,'E_UP',e_up,5.0D0)
3592 call reada(mcmcard,'DELTE',delte,0.1D0)
3593 call readi(mcmcard,'NSWEEP',nsweep,5)
3594 call readi(mcmcard,'NSTEPH',nsteph,0)
3595 call readi(mcmcard,'NSTEPC',nstepc,0)
3596 call reada(mcmcard,'TMIN',tmin,298.0D0)
3597 call reada(mcmcard,'TMAX',tmax,298.0D0)
3598 call readi(mcmcard,'NWINDOW',nwindow,0)
3599 call readi(mcmcard,'PRINT_MC',print_mc,0)
3600 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3601 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3602 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3603 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3604 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3605 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3606 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3607 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3608 if (nwindow.gt.0) then
3609 allocate(winstart(nwindow)) !!el (maxres)
3610 allocate(winend(nwindow)) !!el
3611 allocate(winlen(nwindow)) !!el
3612 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3614 winlen(i)=winend(i)-winstart(i)+1
3617 if (tmax.lt.tmin) tmax=tmin
3618 if (tmax.eq.tmin) then
3622 if (nstepc.gt.0 .and. nsteph.gt.0) then
3623 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3624 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3626 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3627 ! Probabilities of different move types
3628 sumpro_type(0)=0.0D0
3629 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3630 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3631 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3632 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3633 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3634 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3635 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3637 print *,'i',i,' sumprotype',sumpro_type(i)
3638 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3639 print *,'i',i,' sumprotype',sumpro_type(i)
3642 end subroutine mcmread
3643 !-----------------------------------------------------------------------------
3644 subroutine read_minim
3647 ! implicit real*8 (a-h,o-z)
3648 ! include 'DIMENSIONS'
3649 ! include 'COMMON.MINIM'
3650 ! include 'COMMON.IOUNITS'
3651 ! character(len=80) :: ucase
3652 character(len=320) :: minimcard
3654 ! integer :: ntf,ik,iw_pdb
3655 ! real(kind=8) :: var,ene
3657 call card_concat(minimcard,.true.)
3658 call readi(minimcard,'MAXMIN',maxmin,2000)
3659 call readi(minimcard,'MAXFUN',maxfun,5000)
3660 call readi(minimcard,'MINMIN',minmin,maxmin)
3661 call readi(minimcard,'MINFUN',minfun,maxmin)
3662 call reada(minimcard,'TOLF',tolf,1.0D-2)
3663 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3664 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3665 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3666 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3667 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3668 'Options in energy minimization:'
3669 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3670 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3671 'MinMin:',MinMin,' MinFun:',MinFun,&
3672 ' TolF:',TolF,' RTolF:',RTolF
3674 end subroutine read_minim
3675 !-----------------------------------------------------------------------------
3676 subroutine openunits
3678 use MD_data, only: usampl
3681 use control_data, only:out1file
3682 use control, only: getenv_loc
3683 ! implicit real*8 (a-h,o-z)
3684 ! include 'DIMENSIONS'
3687 character(len=16) :: form,nodename
3688 integer :: nodelen,ierror,npos
3690 ! include 'COMMON.SETUP'
3691 ! include 'COMMON.IOUNITS'
3692 ! include 'COMMON.MD'
3693 ! include 'COMMON.CONTROL'
3694 integer :: lenpre,lenpot,lentmp !,ilen
3696 character(len=3) :: out1file_text !,ucase
3697 character(len=3) :: ll
3700 ! integer :: ntf,ik,iw_pdb
3701 ! real(kind=8) :: var,ene
3703 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3704 call getenv_loc("PREFIX",prefix)
3706 call getenv_loc("POT",pot)
3707 call getenv_loc("DIRTMP",tmpdir)
3708 call getenv_loc("CURDIR",curdir)
3709 call getenv_loc("OUT1FILE",out1file_text)
3710 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3711 out1file_text=ucase(out1file_text)
3712 if (out1file_text(1:1).eq."Y") then
3715 out1file=fg_rank.gt.0
3720 if (lentmp.gt.0) then
3721 write (*,'(80(1h!))')
3722 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3723 write (*,'(80(1h!))')
3724 write (*,*)"All output files will be on node /tmp directory."
3726 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3727 if (me.eq.king) then
3728 write (*,*) "The master node is ",nodename
3729 else if (fg_rank.eq.0) then
3730 write (*,*) "I am the CG slave node ",nodename
3732 write (*,*) "I am the FG slave node ",nodename
3735 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3736 lenpre = lentmp+lenpre+1
3738 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3739 ! Get the names and open the input files
3740 #if defined(WINIFL) || defined(WINPGI)
3741 open(1,file=pref_orig(:ilen(pref_orig))// &
3742 '.inp',status='old',readonly,shared)
3743 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3744 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3745 ! Get parameter filenames and open the parameter files.
3746 call getenv_loc('BONDPAR',bondname)
3747 open (ibond,file=bondname,status='old',readonly,shared)
3748 call getenv_loc('THETPAR',thetname)
3749 open (ithep,file=thetname,status='old',readonly,shared)
3750 call getenv_loc('ROTPAR',rotname)
3751 open (irotam,file=rotname,status='old',readonly,shared)
3752 call getenv_loc('TORPAR',torname)
3753 open (itorp,file=torname,status='old',readonly,shared)
3754 call getenv_loc('TORDPAR',tordname)
3755 open (itordp,file=tordname,status='old',readonly,shared)
3756 call getenv_loc('FOURIER',fouriername)
3757 open (ifourier,file=fouriername,status='old',readonly,shared)
3758 call getenv_loc('ELEPAR',elename)
3759 open (ielep,file=elename,status='old',readonly,shared)
3760 call getenv_loc('SIDEPAR',sidename)
3761 open (isidep,file=sidename,status='old',readonly,shared)
3762 #elif (defined CRAY) || (defined AIX)
3763 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3765 ! print *,"Processor",myrank," opened file 1"
3766 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3767 ! print *,"Processor",myrank," opened file 9"
3768 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3769 ! Get parameter filenames and open the parameter files.
3770 call getenv_loc('BONDPAR',bondname)
3771 open (ibond,file=bondname,status='old',action='read')
3772 ! print *,"Processor",myrank," opened file IBOND"
3773 call getenv_loc('THETPAR',thetname)
3774 open (ithep,file=thetname,status='old',action='read')
3775 ! print *,"Processor",myrank," opened file ITHEP"
3776 call getenv_loc('ROTPAR',rotname)
3777 open (irotam,file=rotname,status='old',action='read')
3778 ! print *,"Processor",myrank," opened file IROTAM"
3779 call getenv_loc('TORPAR',torname)
3780 open (itorp,file=torname,status='old',action='read')
3781 ! print *,"Processor",myrank," opened file ITORP"
3782 call getenv_loc('TORDPAR',tordname)
3783 open (itordp,file=tordname,status='old',action='read')
3784 ! print *,"Processor",myrank," opened file ITORDP"
3785 call getenv_loc('SCCORPAR',sccorname)
3786 open (isccor,file=sccorname,status='old',action='read')
3787 ! print *,"Processor",myrank," opened file ISCCOR"
3788 call getenv_loc('FOURIER',fouriername)
3789 open (ifourier,file=fouriername,status='old',action='read')
3790 ! print *,"Processor",myrank," opened file IFOURIER"
3791 call getenv_loc('ELEPAR',elename)
3792 open (ielep,file=elename,status='old',action='read')
3793 ! print *,"Processor",myrank," opened file IELEP"
3794 call getenv_loc('SIDEPAR',sidename)
3795 open (isidep,file=sidename,status='old',action='read')
3796 call getenv_loc('LIPTRANPAR',liptranname)
3797 open (iliptranpar,file=liptranname,status='old',action='read')
3798 call getenv_loc('TUBEPAR',tubename)
3799 open (itube,file=tubename,status='old',action='read')
3801 ! print *,"Processor",myrank," opened file ISIDEP"
3802 ! print *,"Processor",myrank," opened parameter files"
3804 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3805 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3806 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3807 ! Get parameter filenames and open the parameter files.
3808 call getenv_loc('BONDPAR',bondname)
3809 open (ibond,file=bondname,status='old')
3810 call getenv_loc('THETPAR',thetname)
3811 open (ithep,file=thetname,status='old')
3812 call getenv_loc('ROTPAR',rotname)
3813 open (irotam,file=rotname,status='old')
3814 call getenv_loc('TORPAR',torname)
3815 open (itorp,file=torname,status='old')
3816 call getenv_loc('TORDPAR',tordname)
3817 open (itordp,file=tordname,status='old')
3818 call getenv_loc('SCCORPAR',sccorname)
3819 open (isccor,file=sccorname,status='old')
3820 call getenv_loc('FOURIER',fouriername)
3821 open (ifourier,file=fouriername,status='old')
3822 call getenv_loc('ELEPAR',elename)
3823 open (ielep,file=elename,status='old')
3824 call getenv_loc('SIDEPAR',sidename)
3825 open (isidep,file=sidename,status='old')
3826 call getenv_loc('LIPTRANPAR',liptranname)
3827 open (iliptranpar,file=liptranname,status='old')
3828 call getenv_loc('TUBEPAR',tubename)
3829 open (itube,file=tubename,status='old')
3831 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3833 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3834 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3835 ! Get parameter filenames and open the parameter files.
3836 call getenv_loc('BONDPAR',bondname)
3837 open (ibond,file=bondname,status='old',action='read')
3838 call getenv_loc('THETPAR',thetname)
3839 open (ithep,file=thetname,status='old',action='read')
3840 call getenv_loc('ROTPAR',rotname)
3841 open (irotam,file=rotname,status='old',action='read')
3842 call getenv_loc('TORPAR',torname)
3843 open (itorp,file=torname,status='old',action='read')
3844 call getenv_loc('TORDPAR',tordname)
3845 open (itordp,file=tordname,status='old',action='read')
3846 call getenv_loc('SCCORPAR',sccorname)
3847 open (isccor,file=sccorname,status='old',action='read')
3849 call getenv_loc('THETPARPDB',thetname_pdb)
3850 print *,"thetname_pdb ",thetname_pdb
3851 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3852 print *,ithep_pdb," opened"
3854 call getenv_loc('FOURIER',fouriername)
3855 open (ifourier,file=fouriername,status='old',readonly)
3856 call getenv_loc('ELEPAR',elename)
3857 open (ielep,file=elename,status='old',readonly)
3858 call getenv_loc('SIDEPAR',sidename)
3859 open (isidep,file=sidename,status='old',readonly)
3860 call getenv_loc('LIPTRANPAR',liptranname)
3861 open (iliptranpar,file=liptranname,status='old',action='read')
3862 call getenv_loc('TUBEPAR',tubename)
3863 open (itube,file=tubename,status='old',action='read')
3866 call getenv_loc('ROTPARPDB',rotname_pdb)
3867 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3872 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3873 ! Use -DOLDSCP to use hard-coded constants instead.
3875 call getenv_loc('SCPPAR',scpname)
3876 #if defined(WINIFL) || defined(WINPGI)
3877 open (iscpp,file=scpname,status='old',readonly,shared)
3878 #elif (defined CRAY) || (defined AIX)
3879 open (iscpp,file=scpname,status='old',action='read')
3881 open (iscpp,file=scpname,status='old')
3883 open (iscpp,file=scpname,status='old',action='read')
3886 call getenv_loc('PATTERN',patname)
3887 #if defined(WINIFL) || defined(WINPGI)
3888 open (icbase,file=patname,status='old',readonly,shared)
3889 #elif (defined CRAY) || (defined AIX)
3890 open (icbase,file=patname,status='old',action='read')
3892 open (icbase,file=patname,status='old')
3894 open (icbase,file=patname,status='old',action='read')
3897 ! Open output file only for CG processes
3898 ! print *,"Processor",myrank," fg_rank",fg_rank
3899 if (fg_rank.eq.0) then
3901 if (nodes.eq.1) then
3904 npos = dlog10(dfloat(nodes-1))+1
3906 if (npos.lt.3) npos=3
3907 write (liczba,'(i1)') npos
3908 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3910 write (liczba,form) me
3911 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3912 liczba(:ilen(liczba))
3913 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3915 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3917 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3918 liczba(:ilen(liczba))//'.mol2'
3919 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3920 liczba(:ilen(liczba))//'.stat'
3922 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3923 //liczba(:ilen(liczba))//'.stat')
3924 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3927 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3928 liczba(:ilen(liczba))//'.const'
3933 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3934 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3935 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3936 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3937 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3939 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3941 rest2name=prefix(:ilen(prefix))//'.rst'
3943 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3946 #if defined(AIX) || defined(PGI)
3947 if (me.eq.king .or. .not. out1file) &
3948 open(iout,file=outname,status='unknown')
3950 if (fg_rank.gt.0) then
3951 write (liczba,'(i3.3)') myrank/nfgtasks
3952 write (ll,'(bz,i3.3)') fg_rank
3953 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3958 open(igeom,file=intname,status='unknown',position='append')
3959 open(ipdb,file=pdbname,status='unknown')
3960 open(imol2,file=mol2name,status='unknown')
3961 open(istat,file=statname,status='unknown',position='append')
3963 !1out open(iout,file=outname,status='unknown')
3966 if (me.eq.king .or. .not.out1file) &
3967 open(iout,file=outname,status='unknown')
3969 if (fg_rank.gt.0) then
3970 write (liczba,'(i3.3)') myrank/nfgtasks
3971 write (ll,'(bz,i3.3)') fg_rank
3972 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3977 open(igeom,file=intname,status='unknown',access='append')
3978 open(ipdb,file=pdbname,status='unknown')
3979 open(imol2,file=mol2name,status='unknown')
3980 open(istat,file=statname,status='unknown',access='append')
3982 !1out open(iout,file=outname,status='unknown')
3985 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3986 csa_seed=prefix(:lenpre)//'.CSA.seed'
3987 csa_history=prefix(:lenpre)//'.CSA.history'
3988 csa_bank=prefix(:lenpre)//'.CSA.bank'
3989 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3990 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3991 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3992 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3993 csa_int=prefix(:lenpre)//'.int'
3994 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3995 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3996 csa_in=prefix(:lenpre)//'.CSA.in'
3997 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4000 write (iout,'(80(1h-))')
4001 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4002 write (iout,'(80(1h-))')
4003 write (iout,*) "Input file : ",&
4004 pref_orig(:ilen(pref_orig))//'.inp'
4005 write (iout,*) "Output file : ",&
4006 outname(:ilen(outname))
4008 write (iout,*) "Sidechain potential file : ",&
4009 sidename(:ilen(sidename))
4011 write (iout,*) "SCp potential file : ",&
4012 scpname(:ilen(scpname))
4014 write (iout,*) "Electrostatic potential file : ",&
4015 elename(:ilen(elename))
4016 write (iout,*) "Cumulant coefficient file : ",&
4017 fouriername(:ilen(fouriername))
4018 write (iout,*) "Torsional parameter file : ",&
4019 torname(:ilen(torname))
4020 write (iout,*) "Double torsional parameter file : ",&
4021 tordname(:ilen(tordname))
4022 write (iout,*) "SCCOR parameter file : ",&
4023 sccorname(:ilen(sccorname))
4024 write (iout,*) "Bond & inertia constant file : ",&
4025 bondname(:ilen(bondname))
4026 write (iout,*) "Bending parameter file : ",&
4027 thetname(:ilen(thetname))
4028 write (iout,*) "Rotamer parameter file : ",&
4029 rotname(:ilen(rotname))
4032 write (iout,*) "Thetpdb parameter file : ",&
4033 thetname_pdb(:ilen(thetname_pdb))
4036 write (iout,*) "Threading database : ",&
4037 patname(:ilen(patname))
4039 write (iout,*)" DIRTMP : ",&
4041 write (iout,'(80(1h-))')
4044 end subroutine openunits
4045 !-----------------------------------------------------------------------------
4048 use geometry_data, only: nres,dc
4050 ! implicit real*8 (a-h,o-z)
4051 ! include 'DIMENSIONS'
4052 ! include 'COMMON.CHAIN'
4053 ! include 'COMMON.IOUNITS'
4054 ! include 'COMMON.MD'
4057 ! real(kind=8) :: var,ene
4059 open(irest2,file=rest2name,status='unknown')
4060 read(irest2,*) totT,EK,potE,totE,t_bath
4063 ! AL 4/17/17: Now reading d_t(0,:) too
4065 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4068 ! AL 4/17/17: Now reading d_c(0,:) too
4070 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4073 read (irest2,*) iset
4077 end subroutine readrst
4078 !-----------------------------------------------------------------------------
4079 subroutine read_fragments
4083 use control_data, only:out1file
4086 ! implicit real*8 (a-h,o-z)
4087 ! include 'DIMENSIONS'
4091 ! include 'COMMON.SETUP'
4092 ! include 'COMMON.CHAIN'
4093 ! include 'COMMON.IOUNITS'
4094 ! include 'COMMON.MD'
4095 ! include 'COMMON.CONTROL'
4098 ! real(kind=8) :: var,ene
4100 read(inp,*) nset,nfrag,npair,nfrag_back
4102 !el from module energy
4103 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4104 if(.not.allocated(wfrag_back)) then
4105 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4106 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4108 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4109 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4111 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4112 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4115 if(me.eq.king.or..not.out1file) &
4116 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4117 " nfrag_back",nfrag_back
4119 read(inp,*) mset(iset)
4121 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4123 if(me.eq.king.or..not.out1file) &
4124 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4125 ifrag(2,i,iset), qinfrag(i,iset)
4128 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4130 if(me.eq.king.or..not.out1file) &
4131 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4132 ipair(2,i,iset), qinpair(i,iset)
4135 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4136 wfrag_back(3,i,iset),&
4137 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4138 if(me.eq.king.or..not.out1file) &
4139 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4140 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4144 end subroutine read_fragments
4145 !-----------------------------------------------------------------------------
4147 !-----------------------------------------------------------------------------
4151 ! implicit real*8 (a-h,o-z)
4152 ! include 'DIMENSIONS'
4153 ! include 'COMMON.CSA'
4154 ! include 'COMMON.BANK'
4155 ! include 'COMMON.IOUNITS'
4157 ! integer :: ntf,ik,iw_pdb
4158 ! real(kind=8) :: var,ene
4160 open(icsa_in,file=csa_in,status="old",err=100)
4161 read(icsa_in,*) nconf
4162 read(icsa_in,*) jstart,jend
4163 read(icsa_in,*) nstmax
4164 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4165 read(icsa_in,*) nran0,nran1,irr
4166 read(icsa_in,*) nseed
4167 read(icsa_in,*) ntotal,cut1,cut2
4168 read(icsa_in,*) estop
4169 read(icsa_in,*) icmax,irestart
4170 read(icsa_in,*) ntbankm,dele,difcut
4171 read(icsa_in,*) iref,rmscut,pnccut
4172 read(icsa_in,*) ndiff
4179 end subroutine csa_read
4180 !-----------------------------------------------------------------------------
4181 subroutine initial_write
4184 ! implicit real*8 (a-h,o-z)
4185 ! include 'DIMENSIONS'
4186 ! include 'COMMON.CSA'
4187 ! include 'COMMON.BANK'
4188 ! include 'COMMON.IOUNITS'
4190 ! integer :: ntf,ik,iw_pdb
4191 ! real(kind=8) :: var,ene
4193 open(icsa_seed,file=csa_seed,status="unknown")
4194 write(icsa_seed,*) "seed"
4196 #if defined(AIX) || defined(PGI)
4197 open(icsa_history,file=csa_history,status="unknown",&
4200 open(icsa_history,file=csa_history,status="unknown",&
4203 write(icsa_history,*) nconf
4204 write(icsa_history,*) jstart,jend
4205 write(icsa_history,*) nstmax
4206 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4207 write(icsa_history,*) nran0,nran1,irr
4208 write(icsa_history,*) nseed
4209 write(icsa_history,*) ntotal,cut1,cut2
4210 write(icsa_history,*) estop
4211 write(icsa_history,*) icmax,irestart
4212 write(icsa_history,*) ntbankm,dele,difcut
4213 write(icsa_history,*) iref,rmscut,pnccut
4214 write(icsa_history,*) ndiff
4216 write(icsa_history,*)
4219 open(icsa_bank1,file=csa_bank1,status="unknown")
4220 write(icsa_bank1,*) 0
4224 end subroutine initial_write
4225 !-----------------------------------------------------------------------------
4226 subroutine restart_write
4229 ! implicit real*8 (a-h,o-z)
4230 ! include 'DIMENSIONS'
4231 ! include 'COMMON.IOUNITS'
4232 ! include 'COMMON.CSA'
4233 ! include 'COMMON.BANK'
4235 ! integer :: ntf,ik,iw_pdb
4236 ! real(kind=8) :: var,ene
4238 #if defined(AIX) || defined(PGI)
4239 open(icsa_history,file=csa_history,position="append")
4241 open(icsa_history,file=csa_history,access="append")
4243 write(icsa_history,*)
4244 write(icsa_history,*) "This is restart"
4245 write(icsa_history,*)
4246 write(icsa_history,*) nconf
4247 write(icsa_history,*) jstart,jend
4248 write(icsa_history,*) nstmax
4249 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4250 write(icsa_history,*) nran0,nran1,irr
4251 write(icsa_history,*) nseed
4252 write(icsa_history,*) ntotal,cut1,cut2
4253 write(icsa_history,*) estop
4254 write(icsa_history,*) icmax,irestart
4255 write(icsa_history,*) ntbankm,dele,difcut
4256 write(icsa_history,*) iref,rmscut,pnccut
4257 write(icsa_history,*) ndiff
4258 write(icsa_history,*)
4259 write(icsa_history,*) "irestart is: ", irestart
4261 write(icsa_history,*)
4265 end subroutine restart_write
4266 !-----------------------------------------------------------------------------
4268 !-----------------------------------------------------------------------------
4269 subroutine write_pdb(npdb,titelloc,ee)
4271 ! implicit real*8 (a-h,o-z)
4272 ! include 'DIMENSIONS'
4273 ! include 'COMMON.IOUNITS'
4274 character(len=50) :: titelloc1
4275 character*(*) :: titelloc
4276 character(len=3) :: zahl
4277 character(len=5) :: liczba5
4279 integer :: npdb !,ilen
4283 ! real(kind=8) :: var,ene
4287 if (npdb.lt.1000) then
4288 call numstr(npdb,zahl)
4289 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4291 if (npdb.lt.10000) then
4292 write(liczba5,'(i1,i4)') 0,npdb
4294 write(liczba5,'(i5)') npdb
4296 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4298 call pdbout(ee,titelloc1,ipdb)
4301 end subroutine write_pdb
4302 !-----------------------------------------------------------------------------
4304 !-----------------------------------------------------------------------------
4305 subroutine write_thread_summary
4306 ! Thread the sequence through a database of known structures
4307 use control_data, only: refstr
4309 use energy_data, only: n_ene_comp
4311 ! implicit real*8 (a-h,o-z)
4312 ! include 'DIMENSIONS'
4314 use MPI_data !include 'COMMON.INFO'
4317 ! include 'COMMON.CONTROL'
4318 ! include 'COMMON.CHAIN'
4319 ! include 'COMMON.DBASE'
4320 ! include 'COMMON.INTERACT'
4321 ! include 'COMMON.VAR'
4322 ! include 'COMMON.THREAD'
4323 ! include 'COMMON.FFIELD'
4324 ! include 'COMMON.SBRIDGE'
4325 ! include 'COMMON.HEADER'
4326 ! include 'COMMON.NAMES'
4327 ! include 'COMMON.IOUNITS'
4328 ! include 'COMMON.TIME1'
4330 integer,dimension(maxthread) :: ip
4331 real(kind=8),dimension(0:n_ene) :: energia
4333 integer :: i,j,ii,jj,ipj,ik,kk,ist
4334 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4336 write (iout,'(30x,a/)') &
4337 ' *********** Summary threading statistics ************'
4338 write (iout,'(a)') 'Initial energies:'
4339 write (iout,'(a4,2x,a12,14a14,3a8)') &
4340 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4341 'RMSnat','NatCONT','NNCONT','RMS'
4342 ! Energy sort patterns
4347 enet=ener(n_ene-1,ip(i))
4350 if (ener(n_ene-1,ip(j)).lt.enet) then
4352 enet=ener(n_ene-1,ip(j))
4364 ist=nres_base(2,ii)+ipatt(2,i)
4366 energia(i)=ener0(kk,i)
4368 etot=ener0(n_ene_comp+1,i)
4369 rmsnat=ener0(n_ene_comp+2,i)
4370 rms=ener0(n_ene_comp+3,i)
4371 frac=ener0(n_ene_comp+4,i)
4372 frac_nn=ener0(n_ene_comp+5,i)
4375 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4376 i,str_nam(ii),ist+1,&
4377 (energia(print_order(kk)),kk=1,nprint_ene),&
4378 etot,rmsnat,frac,frac_nn,rms
4380 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4381 i,str_nam(ii),ist+1,&
4382 (energia(print_order(kk)),kk=1,nprint_ene),etot
4385 write (iout,'(//a)') 'Final energies:'
4386 write (iout,'(a4,2x,a12,17a14,3a8)') &
4387 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4388 'RMSnat','NatCONT','NNCONT','RMS'
4392 ist=nres_base(2,ii)+ipatt(2,i)
4394 energia(kk)=ener(kk,ik)
4396 etot=ener(n_ene_comp+1,i)
4397 rmsnat=ener(n_ene_comp+2,i)
4398 rms=ener(n_ene_comp+3,i)
4399 frac=ener(n_ene_comp+4,i)
4400 frac_nn=ener(n_ene_comp+5,i)
4401 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4402 i,str_nam(ii),ist+1,&
4403 (energia(print_order(kk)),kk=1,nprint_ene),&
4404 etot,rmsnat,frac,frac_nn,rms
4406 write (iout,'(/a/)') 'IEXAM array:'
4407 write (iout,'(i5)') nexcl
4409 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4411 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4412 'Max. time for threading step ',max_time_for_thread,&
4413 'Average time for threading step: ',ave_time_for_thread
4415 end subroutine write_thread_summary
4416 !-----------------------------------------------------------------------------
4417 subroutine write_stat_thread(ithread,ipattern,ist)
4419 use energy_data, only: n_ene_comp
4421 ! implicit real*8 (a-h,o-z)
4422 ! include "DIMENSIONS"
4423 ! include "COMMON.CONTROL"
4424 ! include "COMMON.IOUNITS"
4425 ! include "COMMON.THREAD"
4426 ! include "COMMON.FFIELD"
4427 ! include "COMMON.DBASE"
4428 ! include "COMMON.NAMES"
4429 real(kind=8),dimension(0:n_ene) :: energia
4431 integer :: ithread,ipattern,ist,i
4432 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4434 #if defined(AIX) || defined(PGI)
4435 open(istat,file=statname,position='append')
4437 open(istat,file=statname,access='append')
4440 energia(i)=ener(i,ithread)
4442 etot=ener(n_ene_comp+1,ithread)
4443 rmsnat=ener(n_ene_comp+2,ithread)
4444 rms=ener(n_ene_comp+3,ithread)
4445 frac=ener(n_ene_comp+4,ithread)
4446 frac_nn=ener(n_ene_comp+5,ithread)
4447 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4448 ithread,str_nam(ipattern),ist+1,&
4449 (energia(print_order(i)),i=1,nprint_ene),&
4450 etot,rmsnat,frac,frac_nn,rms
4453 end subroutine write_stat_thread
4454 !-----------------------------------------------------------------------------
4456 !-----------------------------------------------------------------------------
4457 end module io_config