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,&
747 integer :: ichir1,ichir2
748 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
749 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
750 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 ! For printing parameters after they are read set the following in the UNRES
755 ! setenv PRINT_PARM YES
757 ! To print parameters in LaTeX format rather than as ASCII tables:
761 call getenv_loc("PRINT_PARM",lancuch)
762 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
763 call getenv_loc("LATEX",lancuch)
764 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
766 dwa16=2.0d0**(1.0d0/6.0d0)
768 ! Assign virtual-bond length
771 vblinv2=vblinv*vblinv
773 ! Read the virtual-bond parameters, masses, and moments of inertia
774 ! and Stokes' radii of the peptide group and side chains
776 allocate(dsc(ntyp1)) !(ntyp1)
777 allocate(dsc_inv(ntyp1)) !(ntyp1)
778 allocate(nbondterm(ntyp)) !(ntyp)
779 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
780 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(msc(ntyp+1)) !(ntyp+1)
782 allocate(isc(ntyp+1)) !(ntyp+1)
783 allocate(restok(ntyp+1)) !(ntyp+1)
784 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
790 read (ibond,*) vbldp0,akp,mp,ip,pstok
793 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
794 dsc(i) = vbldsc0(1,i)
798 dsc_inv(i)=1.0D0/dsc(i)
802 read (ibond,*) junk,vbldpDUM,vbldp0,akp,rjunk,mp,ip,pstok
804 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
805 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
806 dsc(i) = vbldsc0(1,i)
810 dsc_inv(i)=1.0D0/dsc(i)
815 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
816 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
818 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
820 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
821 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
823 write (iout,'(13x,3f10.5)') &
824 vbldsc0(j,i),aksc(j,i),abond0(j,i)
828 !----------------------------------------------------
829 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
830 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
831 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
832 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
833 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
834 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
847 ! Read the parameters of the probability distribution/energy expression
848 ! of the virtual-bond valence angles theta
851 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
852 (bthet(j,i,1,1),j=1,2)
853 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
854 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
855 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
859 athet(1,i,1,-1)=athet(1,i,1,1)
860 athet(2,i,1,-1)=athet(2,i,1,1)
861 bthet(1,i,1,-1)=-bthet(1,i,1,1)
862 bthet(2,i,1,-1)=-bthet(2,i,1,1)
863 athet(1,i,-1,1)=-athet(1,i,1,1)
864 athet(2,i,-1,1)=-athet(2,i,1,1)
865 bthet(1,i,-1,1)=bthet(1,i,1,1)
866 bthet(2,i,-1,1)=bthet(2,i,1,1)
870 athet(1,i,-1,-1)=athet(1,-i,1,1)
871 athet(2,i,-1,-1)=-athet(2,-i,1,1)
872 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
873 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
874 athet(1,i,-1,1)=athet(1,-i,1,1)
875 athet(2,i,-1,1)=-athet(2,-i,1,1)
876 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
877 bthet(2,i,-1,1)=bthet(2,-i,1,1)
878 athet(1,i,1,-1)=-athet(1,-i,1,1)
879 athet(2,i,1,-1)=athet(2,-i,1,1)
880 bthet(1,i,1,-1)=bthet(1,-i,1,1)
881 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
886 polthet(j,i)=polthet(j,-i)
889 gthet(j,i)=gthet(j,-i)
897 'Parameters of the virtual-bond valence angles:'
898 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
899 ' ATHETA0 ',' A1 ',' A2 ',&
902 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
903 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
905 write (iout,'(/a/9x,5a/79(1h-))') &
906 'Parameters of the expression for sigma(theta_c):',&
907 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
908 ' ALPH3 ',' SIGMA0C '
910 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
911 (polthet(j,i),j=0,3),sigc0(i)
913 write (iout,'(/a/9x,5a/79(1h-))') &
914 'Parameters of the second gaussian:',&
915 ' THETA0 ',' SIGMA0 ',' G1 ',&
918 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
919 sig0(i),(gthet(j,i),j=1,3)
923 'Parameters of the virtual-bond valence angles:'
924 write (iout,'(/a/9x,5a/79(1h-))') &
925 'Coefficients of expansion',&
926 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
927 ' b1*10^1 ',' b2*10^1 '
929 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
930 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
931 (10*bthet(j,i,1,1),j=1,2)
933 write (iout,'(/a/9x,5a/79(1h-))') &
934 'Parameters of the expression for sigma(theta_c):',&
935 ' alpha0 ',' alph1 ',' alph2 ',&
936 ' alhp3 ',' sigma0c '
938 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
939 (polthet(j,i),j=0,3),sigc0(i)
941 write (iout,'(/a/9x,5a/79(1h-))') &
942 'Parameters of the second gaussian:',&
943 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
946 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
947 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
953 ! Read the parameters of Utheta determined from ab initio surfaces
954 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
956 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
957 ntheterm3,nsingle,ndouble
958 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
960 !----------------------------------------------------
961 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
962 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
963 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
964 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
965 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
966 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
967 !(maxtheterm,-maxthetyp1:maxthetyp1,&
968 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
969 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
970 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
971 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
972 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
973 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
974 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
975 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
976 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
977 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
978 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
979 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
980 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
981 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
982 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
983 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
984 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
986 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
988 ithetyp(i)=-ithetyp(-i)
991 aa0thet(:,:,:,:)=0.0d0
992 aathet(:,:,:,:,:)=0.0d0
993 bbthet(:,:,:,:,:,:)=0.0d0
994 ccthet(:,:,:,:,:,:)=0.0d0
995 ddthet(:,:,:,:,:,:)=0.0d0
996 eethet(:,:,:,:,:,:)=0.0d0
997 ffthet(:,:,:,:,:,:,:)=0.0d0
998 ggthet(:,:,:,:,:,:,:)=0.0d0
1000 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1002 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1003 ! VAR:1=non-glicyne non-proline 2=proline
1004 ! VAR:negative values for D-aminoacid
1006 do j=-nthetyp,nthetyp
1007 do k=-nthetyp,nthetyp
1008 read (ithep,'(6a)',end=111,err=111) res1
1009 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1010 ! VAR: aa0thet is variable describing the average value of Foureir
1011 ! VAR: expansion series
1012 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1013 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1014 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1015 read (ithep,*,end=111,err=111) &
1016 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1017 read (ithep,*,end=111,err=111) &
1018 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1019 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1020 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1021 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1023 read (ithep,*,end=111,err=111) &
1024 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1025 ffthet(lll,llll,ll,i,j,k,iblock),&
1026 ggthet(llll,lll,ll,i,j,k,iblock),&
1027 ggthet(lll,llll,ll,i,j,k,iblock),&
1028 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1033 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1034 ! coefficients of theta-and-gamma-dependent terms are zero.
1035 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1036 ! RECOMENTDED AFTER VERSION 3.3)
1040 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1041 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1043 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1044 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1047 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1049 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1052 ! AND COMMENT THE LOOPS BELOW
1056 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1057 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1059 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1060 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1063 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1065 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1070 ! Substitution for D aminoacids from symmetry.
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1077 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1081 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1082 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1083 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1084 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1090 ffthet(llll,lll,ll,i,j,k,iblock)= &
1091 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1092 ffthet(lll,llll,ll,i,j,k,iblock)= &
1093 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1094 ggthet(llll,lll,ll,i,j,k,iblock)= &
1095 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1096 ggthet(lll,llll,ll,i,j,k,iblock)= &
1097 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1106 ! Control printout of the coefficients of virtual-bond-angle potentials
1109 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1114 write (iout,'(//4a)') &
1115 'Type ',onelett(i),onelett(j),onelett(k)
1116 write (iout,'(//a,10x,a)') " l","a[l]"
1117 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1118 write (iout,'(i2,1pe15.5)') &
1119 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1121 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1122 "b",l,"c",l,"d",l,"e",l
1124 write (iout,'(i2,4(1pe15.5))') m,&
1125 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1126 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1130 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1131 "f+",l,"f-",l,"g+",l,"g-",l
1134 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1135 ffthet(n,m,l,i,j,k,iblock),&
1136 ffthet(m,n,l,i,j,k,iblock),&
1137 ggthet(n,m,l,i,j,k,iblock),&
1138 ggthet(m,n,l,i,j,k,iblock)
1148 write (2,*) "Start reading THETA_PDB",ithep_pdb
1150 ! write (2,*) 'i=',i
1151 read (ithep_pdb,*,err=111,end=111) &
1152 a0thet(i),(athet(j,i,1,1),j=1,2),&
1153 (bthet(j,i,1,1),j=1,2)
1154 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1155 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1156 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1157 sigc0(i)=sigc0(i)**2
1160 athet(1,i,1,-1)=athet(1,i,1,1)
1161 athet(2,i,1,-1)=athet(2,i,1,1)
1162 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1163 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1164 athet(1,i,-1,1)=-athet(1,i,1,1)
1165 athet(2,i,-1,1)=-athet(2,i,1,1)
1166 bthet(1,i,-1,1)=bthet(1,i,1,1)
1167 bthet(2,i,-1,1)=bthet(2,i,1,1)
1170 a0thet(i)=a0thet(-i)
1171 athet(1,i,-1,-1)=athet(1,-i,1,1)
1172 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1173 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1174 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1175 athet(1,i,-1,1)=athet(1,-i,1,1)
1176 athet(2,i,-1,1)=-athet(2,-i,1,1)
1177 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1178 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1179 athet(1,i,1,-1)=-athet(1,-i,1,1)
1180 athet(2,i,1,-1)=athet(2,-i,1,1)
1181 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1182 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1183 theta0(i)=theta0(-i)
1187 polthet(j,i)=polthet(j,-i)
1190 gthet(j,i)=gthet(j,-i)
1193 write (2,*) "End reading THETA_PDB"
1198 !-------------------------------------------
1199 allocate(nlob(ntyp1)) !(ntyp1)
1200 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1201 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1202 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1209 gaussc(:,:,:,:)=0.0D0
1213 ! Read the parameters of the probability distribution/energy expression
1214 ! of the side chains.
1217 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1221 dsc_inv(i)=1.0D0/dsc(i)
1232 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1233 ((blower(k,l,1),l=1,k),k=1,3)
1234 censc(1,1,-i)=censc(1,1,i)
1235 censc(2,1,-i)=censc(2,1,i)
1236 censc(3,1,-i)=-censc(3,1,i)
1238 read (irotam,*,end=112,err=112) bsc(j,i)
1239 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1240 ((blower(k,l,j),l=1,k),k=1,3)
1241 censc(1,j,-i)=censc(1,j,i)
1242 censc(2,j,-i)=censc(2,j,i)
1243 censc(3,j,-i)=-censc(3,j,i)
1244 ! BSC is amplitude of Gaussian
1251 akl=akl+blower(k,m,j)*blower(l,m,j)
1255 if (((k.eq.3).and.(l.ne.3)) &
1256 .or.((l.eq.3).and.(k.ne.3))) then
1257 gaussc(k,l,j,-i)=-akl
1258 gaussc(l,k,j,-i)=-akl
1260 gaussc(k,l,j,-i)=akl
1261 gaussc(l,k,j,-i)=akl
1270 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1273 if (nlobi.gt.0) then
1275 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1276 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1277 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1278 'log h',(bsc(j,i),j=1,nlobi)
1279 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1280 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1282 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1283 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1286 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1287 write (iout,'(a,f10.4,4(16x,f10.4))') &
1288 'Center ',(bsc(j,i),j=1,nlobi)
1289 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1298 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1299 ! added by Urszula Kozlowska 07/11/2007
1301 !el Maximum number of SC local term fitting function coefficiants
1302 !el integer,parameter :: maxsccoef=65
1304 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1307 read (irotam,*,end=112,err=112)
1309 read (irotam,*,end=112,err=112)
1312 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1317 ! Read the parameters of the probability distribution/energy expression
1318 ! of the side chains.
1320 write (2,*) "Start reading ROTAM_PDB"
1322 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1326 dsc_inv(i)=1.0D0/dsc(i)
1337 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1338 ((blower(k,l,1),l=1,k),k=1,3)
1340 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1341 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1342 ((blower(k,l,j),l=1,k),k=1,3)
1349 akl=akl+blower(k,m,j)*blower(l,m,j)
1359 write (2,*) "End reading ROTAM_PDB"
1365 ! Read torsional parameters in old format
1367 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1369 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1370 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1371 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1373 !el from energy module--------
1374 allocate(v1(nterm_old,ntortyp,ntortyp))
1375 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1376 !el---------------------------
1381 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1387 write (iout,'(/a/)') 'Torsional constants:'
1390 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1391 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1397 ! Read torsional parameters
1399 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1400 read (itorp,*,end=113,err=113) ntortyp
1401 !el from energy module---------
1402 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1403 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1405 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1406 allocate(vlor2(maxlor,ntortyp,ntortyp))
1407 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1408 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1410 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1411 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1412 !el---------------------------
1415 !el---------------------------
1417 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1419 itortyp(i)=-itortyp(-i)
1421 ! itortyp(ntyp1)=ntortyp
1422 ! itortyp(-ntyp1)=-ntortyp
1424 write (iout,*) 'ntortyp',ntortyp
1426 do j=-ntortyp+1,ntortyp-1
1427 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1429 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1430 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1433 do k=1,nterm(i,j,iblock)
1434 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1436 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1437 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1438 v0ij=v0ij+si*v1(k,i,j,iblock)
1440 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1441 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1442 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1444 do k=1,nlor(i,j,iblock)
1445 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1446 vlor2(k,i,j),vlor3(k,i,j)
1447 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1450 v0(-i,-j,iblock)=v0ij
1456 write (iout,'(/a/)') 'Torsional constants:'
1458 do i=-ntortyp,ntortyp
1459 do j=-ntortyp,ntortyp
1460 write (iout,*) 'ityp',i,' jtyp',j
1461 write (iout,*) 'Fourier constants'
1462 do k=1,nterm(i,j,iblock)
1463 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1466 write (iout,*) 'Lorenz constants'
1467 do k=1,nlor(i,j,iblock)
1468 write (iout,'(3(1pe15.5))') &
1469 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1475 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1477 ! 6/23/01 Read parameters for double torsionals
1479 !el from energy module------------
1480 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1481 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1482 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1483 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1484 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1485 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1486 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1487 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1488 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1489 !---------------------------------
1493 do j=-ntortyp+1,ntortyp-1
1494 do k=-ntortyp+1,ntortyp-1
1495 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1496 ! write (iout,*) "OK onelett",
1499 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1500 .or. t3.ne.toronelet(k)) then
1501 write (iout,*) "Error in double torsional parameter file",&
1504 call MPI_Finalize(Ierror)
1506 stop "Error in double torsional parameter file"
1508 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1509 ntermd_2(i,j,k,iblock)
1510 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1511 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1512 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1513 ntermd_1(i,j,k,iblock))
1514 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1515 ntermd_1(i,j,k,iblock))
1516 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1517 ntermd_1(i,j,k,iblock))
1518 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1519 ntermd_1(i,j,k,iblock))
1520 ! Martix of D parameters for one dimesional foureir series
1521 do l=1,ntermd_1(i,j,k,iblock)
1522 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1523 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1524 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1525 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1526 ! write(iout,*) "whcodze" ,
1527 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1529 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1530 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1531 v2s(m,l,i,j,k,iblock),&
1532 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1533 ! Martix of D parameters for two dimesional fourier series
1534 do l=1,ntermd_2(i,j,k,iblock)
1536 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1537 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1538 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1539 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1548 write (iout,*) 'Constants for double torsionals'
1551 do j=-ntortyp+1,ntortyp-1
1552 do k=-ntortyp+1,ntortyp-1
1553 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1554 ' nsingle',ntermd_1(i,j,k,iblock),&
1555 ' ndouble',ntermd_2(i,j,k,iblock)
1557 write (iout,*) 'Single angles:'
1558 do l=1,ntermd_1(i,j,k,iblock)
1559 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1560 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1561 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1562 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1565 write (iout,*) 'Pairs of angles:'
1566 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1567 do l=1,ntermd_2(i,j,k,iblock)
1568 write (iout,'(i5,20f10.5)') &
1569 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1572 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1573 do l=1,ntermd_2(i,j,k,iblock)
1574 write (iout,'(i5,20f10.5)') &
1575 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1576 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1585 ! Read of Side-chain backbone correlation parameters
1586 ! Modified 11 May 2012 by Adasko
1589 read (isccor,*,end=119,err=119) nsccortyp
1591 !el from module energy-------------
1592 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1593 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1594 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1595 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1596 !-----------------------------------
1598 !el from module energy-------------
1599 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1601 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1603 isccortyp(i)=-isccortyp(-i)
1605 iscprol=isccortyp(20)
1606 ! write (iout,*) 'ntortyp',ntortyp
1608 !c maxinter is maximum interaction sites
1609 !el from module energy---------
1610 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1611 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1612 -nsccortyp:nsccortyp))
1613 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1614 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1615 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1616 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1617 !-----------------------------------
1621 read (isccor,*,end=119,err=119) &
1622 nterm_sccor(i,j),nlor_sccor(i,j)
1628 nterm_sccor(-i,j)=nterm_sccor(i,j)
1629 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1630 nterm_sccor(i,-j)=nterm_sccor(i,j)
1631 do k=1,nterm_sccor(i,j)
1632 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1634 if (j.eq.iscprol) then
1635 if (i.eq.isccortyp(10)) then
1636 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1637 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1639 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1640 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1641 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1642 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1643 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1644 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1645 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1646 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1649 if (i.eq.isccortyp(10)) then
1650 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1651 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1653 if (j.eq.isccortyp(10)) then
1654 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1655 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1657 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1658 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1659 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1660 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1661 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1662 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1666 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1667 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1668 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1669 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1672 do k=1,nlor_sccor(i,j)
1673 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1674 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1675 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1676 (1+vlor3sccor(k,i,j)**2)
1678 v0sccor(l,i,j)=v0ijsccor
1679 v0sccor(l,-i,j)=v0ijsccor1
1680 v0sccor(l,i,-j)=v0ijsccor2
1681 v0sccor(l,-i,-j)=v0ijsccor3
1687 !el from module energy-------------
1688 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1690 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1691 ! write (iout,*) 'ntortyp',ntortyp
1693 !c maxinter is maximum interaction sites
1694 !el from module energy---------
1695 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1696 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1697 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1698 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1699 !-----------------------------------
1703 read (isccor,*,end=119,err=119) &
1704 nterm_sccor(i,j),nlor_sccor(i,j)
1708 do k=1,nterm_sccor(i,j)
1709 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1711 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1714 do k=1,nlor_sccor(i,j)
1715 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1716 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1717 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1718 (1+vlor3sccor(k,i,j)**2)
1720 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1728 write (iout,'(/a/)') 'Torsional constants:'
1731 write (iout,*) 'ityp',i,' jtyp',j
1732 write (iout,*) 'Fourier constants'
1733 do k=1,nterm_sccor(i,j)
1734 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1736 write (iout,*) 'Lorenz constants'
1737 do k=1,nlor_sccor(i,j)
1738 write (iout,'(3(1pe15.5))') &
1739 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1746 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1747 ! interaction energy of the Gly, Ala, and Pro prototypes.
1752 write (iout,*) "Coefficients of the cumulants"
1754 read (ifourier,*) nloctyp
1755 !write(iout,*) "nloctyp",nloctyp
1756 !el from module energy-------
1757 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1758 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1759 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1760 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1761 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1762 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1763 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1764 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1768 !--------------------------------
1771 read (ifourier,*,end=115,err=115)
1772 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1774 write (iout,*) 'Type',i
1775 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1785 B1tilde(1,-i) =-b(3)
1787 ! b1tilde(1,i)=0.0d0
1788 ! b1tilde(2,i)=0.0d0
1813 Ctilde(1,2,-i)=-b(9)
1817 ! Ctilde(1,1,i)=0.0d0
1818 ! Ctilde(1,2,i)=0.0d0
1819 ! Ctilde(2,1,i)=0.0d0
1820 ! Ctilde(2,2,i)=0.0d0
1838 Dtilde(1,2,-i)=-b(8)
1842 ! Dtilde(1,1,i)=0.0d0
1843 ! Dtilde(1,2,i)=0.0d0
1844 ! Dtilde(2,1,i)=0.0d0
1845 ! Dtilde(2,2,i)=0.0d0
1846 EE(1,1,i)= b(10)+b(11)
1847 EE(2,2,i)=-b(10)+b(11)
1848 EE(2,1,i)= b(12)-b(13)
1849 EE(1,2,i)= b(12)+b(13)
1850 EE(1,1,-i)= b(10)+b(11)
1851 EE(2,2,-i)=-b(10)+b(11)
1852 EE(2,1,-i)=-b(12)+b(13)
1853 EE(1,2,-i)=-b(12)-b(13)
1859 ! ee(2,1,i)=ee(1,2,i)
1863 write (iout,*) 'Type',i
1865 write(iout,*) B1(1,i),B1(2,i)
1867 write(iout,*) B2(1,i),B2(2,i)
1870 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1874 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1878 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1883 ! Read electrostatic-interaction parameters
1888 write (iout,'(/a)') 'Electrostatic interaction constants:'
1889 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1890 'IT','JT','APP','BPP','AEL6','AEL3'
1892 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1893 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1894 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1895 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1900 app (i,j)=epp(i,j)*rri*rri
1901 bpp (i,j)=-2.0D0*epp(i,j)*rri
1902 ael6(i,j)=elpp6(i,j)*4.2D0**6
1903 ael3(i,j)=elpp3(i,j)*4.2D0**3
1905 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1911 ! Read side-chain interaction parameters.
1913 !el from module energy - COMMON.INTERACT-------
1914 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1915 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1916 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1917 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1918 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1927 !--------------------------------
1929 read (isidep,*,end=117,err=117) ipot,expon
1930 if (ipot.lt.1 .or. ipot.gt.5) then
1931 write (iout,'(2a)') 'Error while reading SC interaction',&
1932 'potential file - unknown potential type.'
1934 call MPI_Finalize(Ierror)
1940 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1941 ', exponents are ',expon,2*expon
1942 ! goto (10,20,30,30,40) ipot
1944 !----------------------- LJ potential ---------------------------------
1946 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1947 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1948 (sigma0(i),i=1,ntyp)
1950 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1951 write (iout,'(a/)') 'The epsilon array:'
1952 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1953 write (iout,'(/a)') 'One-body parameters:'
1954 write (iout,'(a,4x,a)') 'residue','sigma'
1955 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1958 !----------------------- LJK potential --------------------------------
1960 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1962 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1964 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1965 write (iout,'(a/)') 'The epsilon array:'
1966 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1967 write (iout,'(/a)') 'One-body parameters:'
1968 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1969 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1973 !---------------------- GB or BP potential -----------------------------
1977 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1979 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1980 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1981 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1982 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1983 ! For the GB potential convert sigma'**2 into chi'
1986 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1990 write (iout,'(/a/)') 'Parameters of the BP potential:'
1991 write (iout,'(a/)') 'The epsilon array:'
1992 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1993 write (iout,'(/a)') 'One-body parameters:'
1994 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1996 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1997 chip(i),alp(i),i=1,ntyp)
2000 !--------------------- GBV potential -----------------------------------
2002 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2003 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2004 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2005 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2007 write (iout,'(/a/)') 'Parameters of the GBV 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,5a)') 'residue',' sigma ',' r0 ',&
2012 's||/s_|_^2',' chip ',' alph '
2013 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2014 sigii(i),chip(i),alp(i),i=1,ntyp)
2017 write(iout,*)"Wrong ipot"
2022 !-----------------------------------------------------------------------
2023 ! Calculate the "working" parameters of SC interactions.
2025 !el from module energy - COMMON.INTERACT-------
2026 allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2027 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2034 !--------------------------------
2043 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2044 sigma(j,i)=sigma(i,j)
2045 rs0(i,j)=dwa16*sigma(i,j)
2049 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2050 'Working parameters of the SC interactions:',&
2051 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2056 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2065 sigeps=dsign(1.0D0,epsij)
2067 aa(i,j)=epsij*rrij*rrij
2068 bb(i,j)=-sigeps*epsij*rrij
2072 sigt1sq=sigma0(i)**2
2073 sigt2sq=sigma0(j)**2
2076 ratsig1=sigt2sq/sigt1sq
2077 ratsig2=1.0D0/ratsig1
2078 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2079 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2080 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2084 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2085 sigmaii(i,j)=rsum_max
2086 sigmaii(j,i)=rsum_max
2088 ! sigmaii(i,j)=r0(i,j)
2089 ! sigmaii(j,i)=r0(i,j)
2091 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2092 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2093 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2094 augm(i,j)=epsij*r_augm**(2*expon)
2095 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2102 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2103 restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
2104 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2109 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2114 ! Define the SC-p interaction constants (hard-coded; old style)
2117 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2119 ! aad(i,1)=0.3D0*4.0D0**12
2120 ! Following line for constants currently implemented
2121 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2122 aad(i,1)=1.5D0*4.0D0**12
2123 ! aad(i,1)=0.17D0*5.6D0**12
2125 ! "Soft" SC-p repulsion
2127 ! Following line for constants currently implemented
2128 ! aad(i,1)=0.3D0*4.0D0**6
2129 ! "Hard" SC-p repulsion
2130 bad(i,1)=3.0D0*4.0D0**6
2131 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2140 ! 8/9/01 Read the SC-p interaction constants from file
2143 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2146 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2147 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2148 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2149 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2153 write (iout,*) "Parameters of SC-p interactions:"
2155 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2156 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2162 ! Define the constants of the disulfide bridge
2166 ! Old arbitrary potential - commented out.
2171 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2172 ! energy surface of diethyl disulfide.
2173 ! A. Liwo and U. Kozlowska, 11/24/03
2190 write (iout,'(/a)') "Disulfide bridge parameters:"
2191 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2192 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2193 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2194 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2198 111 write (iout,*) "Error reading bending energy parameters."
2200 112 write (iout,*) "Error reading rotamer energy parameters."
2202 113 write (iout,*) "Error reading torsional energy parameters."
2204 114 write (iout,*) "Error reading double torsional energy parameters."
2206 115 write (iout,*) &
2207 "Error reading cumulant (multibody energy) parameters."
2209 116 write (iout,*) "Error reading electrostatic energy parameters."
2211 117 write (iout,*) "Error reading side chain interaction parameters."
2213 118 write (iout,*) "Error reading SCp interaction parameters."
2215 119 write (iout,*) "Error reading SCCOR parameters"
2218 call MPI_Finalize(Ierror)
2222 end subroutine parmread
2224 !-----------------------------------------------------------------------------
2226 !-----------------------------------------------------------------------------
2227 subroutine printmat(ldim,m,n,iout,key,a)
2230 character(len=3),dimension(n) :: key
2231 real(kind=8),dimension(ldim,n) :: a
2233 integer :: i,j,k,m,iout,nlim
2237 write (iout,1000) (key(k),k=i,nlim)
2239 1000 format (/5x,8(6x,a3))
2240 1020 format (/80(1h-)/)
2242 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2245 1010 format (a3,2x,8(f9.4))
2247 end subroutine printmat
2248 !-----------------------------------------------------------------------------
2250 !-----------------------------------------------------------------------------
2252 ! Read the PDB file and convert the peptide geometry into virtual-chain
2255 use energy_data, only: itype
2259 use control, only: rescode
2260 ! implicit real*8 (a-h,o-z)
2261 ! include 'DIMENSIONS'
2262 ! include 'COMMON.LOCAL'
2263 ! include 'COMMON.VAR'
2264 ! include 'COMMON.CHAIN'
2265 ! include 'COMMON.INTERACT'
2266 ! include 'COMMON.IOUNITS'
2267 ! include 'COMMON.GEO'
2268 ! include 'COMMON.NAMES'
2269 ! include 'COMMON.CONTROL'
2270 ! include 'COMMON.DISTFIT'
2271 ! include 'COMMON.SETUP'
2272 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2274 logical :: lprn=.true.,fail
2275 real(kind=8),dimension(3) :: e1,e2,e3
2276 real(kind=8) :: dcj,efree_temp
2277 character(len=3) :: seq,res
2278 character(len=5) :: atom
2279 character(len=80) :: card
2280 real(kind=8),dimension(3,20) :: sccor
2281 integer :: kkk,lll,icha,kupa !rescode,
2284 integer,dimension(2,maxres/3) :: hfrag_alloc
2285 integer,dimension(4,maxres/3) :: bfrag_alloc
2286 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2292 ! write (2,*) "UNRES_PDB",unres_pdb
2300 !-----------------------------
2301 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2302 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2305 read (ipdbin,'(a80)',end=10) card
2306 ! write (iout,'(a)') card
2307 if (card(:5).eq.'HELIX') then
2310 read(card(22:25),*) hfrag(1,nhfrag)
2311 read(card(34:37),*) hfrag(2,nhfrag)
2313 if (card(:5).eq.'SHEET') then
2316 read(card(24:26),*) bfrag(1,nbfrag)
2317 read(card(35:37),*) bfrag(2,nbfrag)
2318 !rc----------------------------------------
2319 !rc to be corrected !!!
2320 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2321 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2322 !rc----------------------------------------
2324 if (card(:3).eq.'END') then
2326 else if (card(:3).eq.'TER') then
2330 itype(ires_old)=ntyp1
2331 itype(ires_old-1)=ntyp1
2333 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2336 dc(j,ires)=sccor(j,iii)
2339 call sccenter(ires,iii,sccor)
2344 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2345 ! Fish out the ATOM cards.
2346 if (index(card(1:4),'ATOM').gt.0) then
2347 read (card(12:16),*) atom
2348 ! write (iout,*) "! ",atom," !",ires
2349 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2350 read (card(23:26),*) ires
2351 read (card(18:20),'(a3)') res
2352 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2353 ! & " ires_old",ires_old
2354 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2355 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2356 if (ires-ishift+ishift1.ne.ires_old) then
2357 ! Calculate the CM of the preceding residue.
2358 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2360 ! write (iout,*) "Calculating sidechain center iii",iii
2363 dc(j,ires+nres)=sccor(j,iii)
2366 call sccenter(ires_old,iii,sccor)
2370 ! Start new residue.
2371 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2374 else if (ibeg.eq.1) then
2375 write (iout,*) "BEG ires",ires
2377 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2381 ires=ires-ishift+ishift1
2383 ! write (iout,*) "ishift",ishift," ires",ires,&
2384 ! " ires_old",ires_old
2386 else if (ibeg.eq.2) then
2388 ishift=-ires_old+ires-1 !!!!!
2389 ishift1=ishift1-1 !!!!!
2390 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2391 ires=ires-ishift+ishift1
2395 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2396 ires=ires-ishift+ishift1
2399 if (res.eq.'ACE' .or. res.eq.'NHE') then
2402 itype(ires)=rescode(ires,res,0)
2405 ires=ires-ishift+ishift1
2407 ! write (iout,*) "ires_old",ires_old," ires",ires
2408 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2411 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2412 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2413 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2414 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2415 ! write (iout,*) "backbone ",atom
2417 write (iout,'(2i3,2x,a,3f8.3)') &
2418 ires,itype(ires),res,(c(j,ires),j=1,3)
2422 sccor(j,iii)=c(j,ires)
2424 ! write (*,*) card(23:27),ires,itype(ires)
2425 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2426 atom.ne.'N' .and. atom.ne.'C' .and. &
2427 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2428 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2429 ! write (iout,*) "sidechain ",atom
2431 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2435 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2436 if (ires.eq.0) return
2437 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2441 ! write (iout,*) i,itype(i)
2442 ! if (itype(i).eq.ntyp1) then
2443 ! write (iout,*) "dummy",i,itype(i)
2445 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2446 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2450 if (itype(i).eq.ntyp1) then
2451 if (itype(i+1).eq.ntyp1) then
2452 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2453 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
2454 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
2456 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2457 ! print *,i,'tu dochodze'
2458 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2466 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2470 dcj=(c(j,i-2)-c(j,i-3))/2.0
2471 if (dcj.eq.0) dcj=1.23591524223
2476 else !itype(i+1).eq.ntyp1
2478 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2479 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2486 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2490 dcj=(c(j,i+3)-c(j,i+2))/2.0
2491 if (dcj.eq.0) dcj=1.23591524223
2496 endif !itype(i+1).eq.ntyp1
2497 endif !itype.eq.ntyp1
2500 ! Calculate the CM of the last side chain.
2504 dc(j,ires)=sccor(j,iii)
2507 call sccenter(ires,iii,sccor)
2513 if (itype(nres).ne.10) then
2517 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2518 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2525 c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
2529 dcj=c(j,nres-2)-c(j,nres-3)
2530 c(j,nres)=c(j,nres-1)+dcj
2531 c(j,2*nres)=c(j,nres)
2535 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2537 if (nres.ne.nres0) then
2538 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2540 stop "Error nres value in WHAM input"
2543 !---------------------------------
2544 !el reallocate tables
2547 ! hfrag_alloc(j,i)=hfrag(j,i)
2550 ! bfrag_alloc(j,i)=bfrag(j,i)
2556 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2557 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2558 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2559 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2563 ! hfrag(j,i)=hfrag_alloc(j,i)
2568 ! bfrag(j,i)=bfrag_alloc(j,i)
2571 !el end reallocate tables
2572 !---------------------------------
2580 c(j,2*nres)=c(j,nres)
2582 if (itype(1).eq.ntyp1) then
2586 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2587 call refsys(2,3,4,e1,e2,e3,fail)
2594 c(j,1)=c(j,2)-3.8d0*e2(j)
2604 ! Copy the coordinates to reference coordinates
2610 ! Calculate internal coordinates.
2612 write (iout,'(/a)') &
2613 "Cartesian coordinates of the reference structure"
2614 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2615 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2617 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2618 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2619 (c(j,ires+nres),j=1,3)
2622 ! znamy już nres wiec mozna alokowac tablice
2623 ! Calculate internal coordinates.
2624 if(me.eq.king.or..not.out1file)then
2625 write (iout,'(a)') &
2626 "Backbone and SC coordinates as read from the PDB"
2628 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2629 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2630 (c(j,nres+ires),j=1,3)
2634 if(.not.allocated(vbld)) then
2635 allocate(vbld(2*nres))
2640 if(.not.allocated(vbld_inv)) then
2641 allocate(vbld_inv(2*nres))
2647 if(.not.allocated(theta)) then
2648 allocate(theta(nres+2))
2652 if(.not.allocated(phi)) allocate(phi(nres+2))
2653 if(.not.allocated(alph)) allocate(alph(nres+2))
2654 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2655 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2656 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2657 if(.not.allocated(costtab)) allocate(costtab(nres))
2658 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2659 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2660 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2661 if(.not.allocated(xxref)) allocate(xxref(nres))
2662 if(.not.allocated(yyref)) allocate(yyref(nres))
2663 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2664 if(.not.allocated(dc_norm)) then
2665 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2666 allocate(dc_norm(3,0:2*nres+2))
2670 call int_from_cart(.true.,.false.)
2671 call sc_loc_geom(.false.)
2673 thetaref(i)=theta(i)
2683 dc(j,i)=c(j,i+1)-c(j,i)
2684 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2689 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2690 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2692 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2696 ! Copy the coordinates to reference coordinates
2697 ! Splits to single chain if occurs
2701 ! cref(j,i,cou)=c(j,i)
2705 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2706 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2707 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2708 !-----------------------------
2714 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2716 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2719 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2724 cref(j,i,cou)=c(j,i)
2725 cref(j,i+nres,cou)=c(j,i+nres)
2727 chain_rep(j,lll,kkk)=c(j,i)
2728 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2732 write (iout,*) chain_length
2733 if (chain_length.eq.0) chain_length=nres
2735 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2736 chain_rep(j,chain_length+nres,symetr) &
2737 =chain_rep(j,chain_length+nres,1)
2740 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2742 ! do kkk=1,chain_length
2743 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2747 ! makes copy of chains
2748 write (iout,*) "symetr", symetr
2750 if (symetr.gt.1) then
2757 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2763 ! write (iout,*) i,icha
2764 do lll=1,chain_length
2766 if (cou.le.nres) then
2768 kupa=mod(lll,chain_length)
2769 iprzes=(kkk-1)*chain_length+lll
2770 if (kupa.eq.0) kupa=chain_length
2771 ! write (iout,*) "kupa", kupa
2772 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2773 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2780 !-koniec robienia kopii
2783 write (iout,*) "nowa struktura", nperm
2785 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2787 cref(3,i,kkk),cref(1,nres+i,kkk),&
2788 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2790 100 format (//' alpha-carbon coordinates ',&
2791 ' centroid coordinates'/ &
2792 ' ', 6X,'X',11X,'Y',11X,'Z', &
2793 10X,'X',11X,'Y',11X,'Z')
2794 110 format (a,'(',i3,')',6f12.5)
2800 bfrag(i,j)=bfrag(i,j)-ishift
2806 hfrag(i,j)=hfrag(i,j)-ishift
2812 end subroutine readpdb
2813 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2814 !-----------------------------------------------------------------------------
2816 !-----------------------------------------------------------------------------
2817 subroutine read_control
2831 use random, only: random_init
2832 ! implicit real*8 (a-h,o-z)
2833 ! include 'DIMENSIONS'
2835 use prng, only:prng_restart
2837 logical :: OKRandom!, prng_restart
2840 ! include 'COMMON.IOUNITS'
2841 ! include 'COMMON.TIME1'
2842 ! include 'COMMON.THREAD'
2843 ! include 'COMMON.SBRIDGE'
2844 ! include 'COMMON.CONTROL'
2845 ! include 'COMMON.MCM'
2846 ! include 'COMMON.MAP'
2847 ! include 'COMMON.HEADER'
2848 ! include 'COMMON.CSA'
2849 ! include 'COMMON.CHAIN'
2850 ! include 'COMMON.MUCA'
2851 ! include 'COMMON.MD'
2852 ! include 'COMMON.FFIELD'
2853 ! include 'COMMON.INTERACT'
2854 ! include 'COMMON.SETUP'
2855 !el integer :: KDIAG,ICORFL,IXDR
2856 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2857 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2858 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2859 ! character(len=80) :: ucase
2860 character(len=640) :: controlcard
2862 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2868 read (INP,'(a)') titel
2869 call card_concat(controlcard,.true.)
2870 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2871 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2872 call reada(controlcard,'SEED',seed,0.0D0)
2873 call random_init(seed)
2874 ! Set up the time limit (caution! The time must be input in minutes!)
2875 read_cart=index(controlcard,'READ_CART').gt.0
2876 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2877 call readi(controlcard,'SYM',symetr,1)
2878 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2879 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2880 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2881 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2882 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2883 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2884 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2885 call reada(controlcard,'DRMS',drms,0.1D0)
2886 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2887 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2888 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2889 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2890 write (iout,'(a,f10.1)')'DRMS = ',drms
2891 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2892 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2894 call readi(controlcard,'NZ_START',nz_start,0)
2895 call readi(controlcard,'NZ_END',nz_end,0)
2896 call readi(controlcard,'IZ_SC',iz_sc,0)
2897 timlim=60.0D0*timlim
2898 safety = 60.0d0*safety
2901 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2902 !C Varibles set size of box
2903 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2904 call reada(controlcard,'BOXY',boxysize,100.0d0)
2905 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2906 ! CUTOFFF ON ELECTROSTATICS
2907 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2908 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2910 !C-------------------------
2911 minim=(index(controlcard,'MINIMIZE').gt.0)
2912 dccart=(index(controlcard,'CART').gt.0)
2913 overlapsc=(index(controlcard,'OVERLAP').gt.0)
2914 overlapsc=.not.overlapsc
2915 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2916 searchsc=.not.searchsc
2917 sideadd=(index(controlcard,'SIDEADD').gt.0)
2918 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2919 outpdb=(index(controlcard,'PDBOUT').gt.0)
2920 outmol2=(index(controlcard,'MOL2OUT').gt.0)
2921 pdbref=(index(controlcard,'PDBREF').gt.0)
2922 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2923 indpdb=index(controlcard,'PDBSTART')
2924 extconf=(index(controlcard,'EXTCONF').gt.0)
2925 call readi(controlcard,'IPRINT',iprint,0)
2926 call readi(controlcard,'MAXGEN',maxgen,10000)
2927 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2928 call readi(controlcard,"KDIAG",kdiag,0)
2929 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2930 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2931 write (iout,*) "RESCALE_MODE",rescale_mode
2932 split_ene=index(controlcard,'SPLIT_ENE').gt.0
2933 if (index(controlcard,'REGULAR').gt.0.0D0) then
2934 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2938 if (index(controlcard,'CHECKGRAD').gt.0) then
2940 if (index(controlcard,'CART').gt.0) then
2942 elseif (index(controlcard,'CARINT').gt.0) then
2947 elseif (index(controlcard,'THREAD').gt.0) then
2949 call readi(controlcard,'THREAD',nthread,0)
2950 if (nthread.gt.0) then
2951 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2954 write (iout,'(a)')'A number has to follow the THREAD keyword.'
2955 stop 'Error termination in Read_Control.'
2957 else if (index(controlcard,'MCMA').gt.0) then
2959 else if (index(controlcard,'MCEE').gt.0) then
2961 else if (index(controlcard,'MULTCONF').gt.0) then
2963 else if (index(controlcard,'MAP').gt.0) then
2965 call readi(controlcard,'MAP',nmap,0)
2966 else if (index(controlcard,'CSA').gt.0) then
2968 !rc else if (index(controlcard,'ZSCORE').gt.0) then
2970 !rc ZSCORE is rm from UNRES, modecalc=9 is available
2973 !fcm else if (index(controlcard,'MCMF').gt.0) then
2975 else if (index(controlcard,'SOFTREG').gt.0) then
2977 else if (index(controlcard,'CHECK_BOND').gt.0) then
2979 else if (index(controlcard,'TEST').gt.0) then
2981 else if (index(controlcard,'MD').gt.0) then
2983 else if (index(controlcard,'RE ').gt.0) then
2987 lmuca=index(controlcard,'MUCA').gt.0
2988 call readi(controlcard,'MUCADYN',mucadyn,0)
2989 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
2990 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
2992 write (iout,*) 'MUCADYN=',mucadyn
2993 write (iout,*) 'MUCASMOOTH=',muca_smooth
2996 iscode=index(controlcard,'ONE_LETTER')
2997 indphi=index(controlcard,'PHI')
2998 indback=index(controlcard,'BACK')
2999 iranconf=index(controlcard,'RAND_CONF')
3000 i2ndstr=index(controlcard,'USE_SEC_PRED')
3001 gradout=index(controlcard,'GRADOUT').gt.0
3002 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3003 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3004 if (me.eq.king .or. .not.out1file ) &
3005 write (iout,*) "DISTCHAINMAX",distchainmax
3007 if(me.eq.king.or..not.out1file) &
3008 write (iout,'(2a)') diagmeth(kdiag),&
3009 ' routine used to diagonalize matrices.'
3011 end subroutine read_control
3012 !-----------------------------------------------------------------------------
3013 subroutine read_REMDpar
3015 ! Read REMD settings
3022 use control_data, only:out1file
3023 ! implicit real*8 (a-h,o-z)
3024 ! include 'DIMENSIONS'
3025 ! include 'COMMON.IOUNITS'
3026 ! include 'COMMON.TIME1'
3027 ! include 'COMMON.MD'
3030 !el include 'COMMON.LANGEVIN'
3032 !el include 'COMMON.LANGEVIN.lang0'
3034 ! include 'COMMON.INTERACT'
3035 ! include 'COMMON.NAMES'
3036 ! include 'COMMON.GEO'
3037 ! include 'COMMON.REMD'
3038 ! include 'COMMON.CONTROL'
3039 ! include 'COMMON.SETUP'
3040 ! character(len=80) :: ucase
3041 character(len=320) :: controlcard
3042 character(len=3200) :: controlcard1
3043 integer :: iremd_m_total
3046 ! real(kind=8) :: var,ene
3048 if(me.eq.king.or..not.out1file) &
3049 write (iout,*) "REMD setup"
3051 call card_concat(controlcard,.true.)
3052 call readi(controlcard,"NREP",nrep,3)
3053 call readi(controlcard,"NSTEX",nstex,1000)
3054 call reada(controlcard,"RETMIN",retmin,10.0d0)
3055 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3056 mremdsync=(index(controlcard,'SYNC').gt.0)
3057 call readi(controlcard,"NSYN",i_sync_step,100)
3058 restart1file=(index(controlcard,'REST1FILE').gt.0)
3059 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3060 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3061 if(max_cache_traj_use.gt.max_cache_traj) &
3062 max_cache_traj_use=max_cache_traj
3063 if(me.eq.king.or..not.out1file) then
3064 !d if (traj1file) then
3065 !rc caching is in testing - NTWX is not ignored
3066 !d write (iout,*) "NTWX value is ignored"
3067 !d write (iout,*) " trajectory is stored to one file by master"
3068 !d write (iout,*) " before exchange at NSTEX intervals"
3070 write (iout,*) "NREP= ",nrep
3071 write (iout,*) "NSTEX= ",nstex
3072 write (iout,*) "SYNC= ",mremdsync
3073 write (iout,*) "NSYN= ",i_sync_step
3074 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3077 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3078 if (index(controlcard,'TLIST').gt.0) then
3080 call card_concat(controlcard1,.true.)
3081 read(controlcard1,*) (remd_t(i),i=1,nrep)
3082 if(me.eq.king.or..not.out1file) &
3083 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3086 if (index(controlcard,'MLIST').gt.0) then
3088 call card_concat(controlcard1,.true.)
3089 read(controlcard1,*) (remd_m(i),i=1,nrep)
3090 if(me.eq.king.or..not.out1file) then
3091 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3094 iremd_m_total=iremd_m_total+remd_m(i)
3096 write (iout,*) 'Total number of replicas ',iremd_m_total
3099 if(me.eq.king.or..not.out1file) &
3100 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3102 end subroutine read_REMDpar
3103 !-----------------------------------------------------------------------------
3104 subroutine read_MDpar
3108 use control_data, only: r_cut,rlamb,out1file
3110 use geometry_data, only: pi
3112 ! implicit real*8 (a-h,o-z)
3113 ! include 'DIMENSIONS'
3114 ! include 'COMMON.IOUNITS'
3115 ! include 'COMMON.TIME1'
3116 ! include 'COMMON.MD'
3119 !el include 'COMMON.LANGEVIN'
3121 !el include 'COMMON.LANGEVIN.lang0'
3123 ! include 'COMMON.INTERACT'
3124 ! include 'COMMON.NAMES'
3125 ! include 'COMMON.GEO'
3126 ! include 'COMMON.SETUP'
3127 ! include 'COMMON.CONTROL'
3128 ! include 'COMMON.SPLITELE'
3129 ! character(len=80) :: ucase
3130 character(len=320) :: controlcard
3135 call card_concat(controlcard,.true.)
3136 call readi(controlcard,"NSTEP",n_timestep,1000000)
3137 call readi(controlcard,"NTWE",ntwe,100)
3138 call readi(controlcard,"NTWX",ntwx,1000)
3139 call reada(controlcard,"DT",d_time,1.0d-1)
3140 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3141 call reada(controlcard,"DAMAX",damax,1.0d1)
3142 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3143 call readi(controlcard,"LANG",lang,0)
3144 RESPA = index(controlcard,"RESPA") .gt. 0
3145 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3146 ntime_split0=ntime_split
3147 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3148 ntime_split0=ntime_split
3149 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3150 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3151 rest = index(controlcard,"REST").gt.0
3152 tbf = index(controlcard,"TBF").gt.0
3153 usampl = index(controlcard,"USAMPL").gt.0
3154 mdpdb = index(controlcard,"MDPDB").gt.0
3155 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3156 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3157 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3158 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3159 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3160 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3161 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3162 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3163 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3164 large = index(controlcard,"LARGE").gt.0
3165 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3166 rattle = index(controlcard,"RATTLE").gt.0
3167 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3173 if(me.eq.king.or..not.out1file) then
3175 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3177 write (iout,'(a)') "The units are:"
3178 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3179 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3180 " acceleration: angstrom/(48.9 fs)**2"
3181 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3183 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3184 write (iout,'(a60,f10.5,a)') &
3185 "Initial time step of numerical integration:",d_time,&
3187 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3189 write (iout,'(2a,i4,a)') &
3190 "A-MTS algorithm used; initial time step for fast-varying",&
3191 " short-range forces split into",ntime_split," steps."
3192 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3193 r_cut," lambda",rlamb
3195 write (iout,'(2a,f10.5)') &
3196 "Maximum acceleration threshold to reduce the time step",&
3197 "/increase split number:",damax
3198 write (iout,'(2a,f10.5)') &
3199 "Maximum predicted energy drift to reduce the timestep",&
3200 "/increase split number:",edriftmax
3201 write (iout,'(a60,f10.5)') &
3202 "Maximum velocity threshold to reduce velocities:",dvmax
3203 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3204 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3205 if (rattle) write (iout,'(a60)') &
3206 "Rattle algorithm used to constrain the virtual bonds"
3210 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3211 call reada(controlcard,"RWAT",rwat,1.4d0)
3212 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3213 surfarea=index(controlcard,"SURFAREA").gt.0
3214 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3215 if(me.eq.king.or..not.out1file)then
3216 write (iout,'(/a,$)') "Langevin dynamics calculation"
3218 write (iout,'(a/)') &
3219 " with direct integration of Langevin equations"
3220 else if (lang.eq.2) then
3221 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3222 else if (lang.eq.3) then
3223 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3224 else if (lang.eq.4) then
3225 write (iout,'(a/)') " in overdamped mode"
3227 write (iout,'(//a,i5)') &
3228 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3231 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3232 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3233 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3234 write (iout,'(a60,f10.5)') &
3235 "Scaling factor of the friction forces:",scal_fric
3236 if (surfarea) write (iout,'(2a,i10,a)') &
3237 "Friction coefficients will be scaled by solvent-accessible",&
3238 " surface area every",reset_fricmat," steps."
3240 ! Calculate friction coefficients and bounds of stochastic forces
3241 eta=6*pi*cPoise*etawat
3242 if(me.eq.king.or..not.out1file) &
3243 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3245 gamp=scal_fric*(pstok+rwat)*eta
3246 stdfp=dsqrt(2*Rb*t_bath/d_time)
3247 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3249 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3250 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3252 if(me.eq.king.or..not.out1file)then
3253 write (iout,'(/2a/)') &
3254 "Radii of site types and friction coefficients and std's of",&
3255 " stochastic forces of fully exposed sites"
3256 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3258 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3259 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3263 if(me.eq.king.or..not.out1file)then
3264 write (iout,'(a)') "Berendsen bath calculation"
3265 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3266 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3268 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3269 count_reset_moment," steps"
3271 write (iout,'(a,i10,a)') &
3272 "Velocities will be reset at random every",count_reset_vel,&
3276 if(me.eq.king.or..not.out1file) &
3277 write (iout,'(a31)') "Microcanonical mode calculation"
3279 if(me.eq.king.or..not.out1file)then
3280 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3282 write(iout,*) "MD running with constraints."
3283 write(iout,*) "Equilibration time ", eq_time, " mtus."
3284 write(iout,*) "Constraining ", nfrag," fragments."
3285 write(iout,*) "Length of each fragment, weight and q0:"
3287 write (iout,*) "Set of restraints #",iset
3289 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3290 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3292 write(iout,*) "constraints between ", npair, "fragments."
3293 write(iout,*) "constraint pairs, weights and q0:"
3295 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3296 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3298 write(iout,*) "angle constraints within ", nfrag_back,&
3299 "backbone fragments."
3300 write(iout,*) "fragment, weights:"
3302 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3303 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3304 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3307 iset=mod(kolor,nset)+1
3310 if(me.eq.king.or..not.out1file) &
3311 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3313 end subroutine read_MDpar
3314 !-----------------------------------------------------------------------------
3318 ! implicit real*8 (a-h,o-z)
3319 ! include 'DIMENSIONS'
3320 ! include 'COMMON.MAP'
3321 ! include 'COMMON.IOUNITS'
3322 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3323 character(len=80) :: mapcard !,ucase
3326 ! real(kind=8) :: var,ene
3329 read (inp,'(a)') mapcard
3330 mapcard=ucase(mapcard)
3331 if (index(mapcard,'PHI').gt.0) then
3333 else if (index(mapcard,'THE').gt.0) then
3335 else if (index(mapcard,'ALP').gt.0) then
3337 else if (index(mapcard,'OME').gt.0) then
3340 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3341 stop 'Error - illegal variable spec in MAP card.'
3343 call readi (mapcard,'RES1',res1(imap),0)
3344 call readi (mapcard,'RES2',res2(imap),0)
3345 if (res1(imap).eq.0) then
3346 res1(imap)=res2(imap)
3347 else if (res2(imap).eq.0) then
3348 res2(imap)=res1(imap)
3350 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3351 write (iout,'(a)') &
3352 'Error - illegal definition of variable group in MAP.'
3353 stop 'Error - illegal definition of variable group in MAP.'
3355 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3356 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3357 call readi(mapcard,'NSTEP',nstep(imap),0)
3358 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3359 write (iout,'(a)') &
3360 'Illegal boundary and/or step size specification in MAP.'
3361 stop 'Illegal boundary and/or step size specification in MAP.'
3365 end subroutine map_read
3366 !-----------------------------------------------------------------------------
3369 use control_data, only: vdisulf
3371 ! implicit real*8 (a-h,o-z)
3372 ! include 'DIMENSIONS'
3373 ! include 'COMMON.IOUNITS'
3374 ! include 'COMMON.GEO'
3375 ! include 'COMMON.CSA'
3376 ! include 'COMMON.BANK'
3377 ! include 'COMMON.CONTROL'
3378 ! character(len=80) :: ucase
3379 character(len=620) :: mcmcard
3381 ! integer :: ntf,ik,iw_pdb
3382 ! real(kind=8) :: var,ene
3384 call card_concat(mcmcard,.true.)
3386 call readi(mcmcard,'NCONF',nconf,50)
3387 call readi(mcmcard,'NADD',nadd,0)
3388 call readi(mcmcard,'JSTART',jstart,1)
3389 call readi(mcmcard,'JEND',jend,1)
3390 call readi(mcmcard,'NSTMAX',nstmax,500000)
3391 call readi(mcmcard,'N0',n0,1)
3392 call readi(mcmcard,'N1',n1,6)
3393 call readi(mcmcard,'N2',n2,4)
3394 call readi(mcmcard,'N3',n3,0)
3395 call readi(mcmcard,'N4',n4,0)
3396 call readi(mcmcard,'N5',n5,0)
3397 call readi(mcmcard,'N6',n6,10)
3398 call readi(mcmcard,'N7',n7,0)
3399 call readi(mcmcard,'N8',n8,0)
3400 call readi(mcmcard,'N9',n9,0)
3401 call readi(mcmcard,'N14',n14,0)
3402 call readi(mcmcard,'N15',n15,0)
3403 call readi(mcmcard,'N16',n16,0)
3404 call readi(mcmcard,'N17',n17,0)
3405 call readi(mcmcard,'N18',n18,0)
3407 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3409 call readi(mcmcard,'NDIFF',ndiff,2)
3410 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3411 call readi(mcmcard,'IS1',is1,1)
3412 call readi(mcmcard,'IS2',is2,8)
3413 call readi(mcmcard,'NRAN0',nran0,4)
3414 call readi(mcmcard,'NRAN1',nran1,2)
3415 call readi(mcmcard,'IRR',irr,1)
3416 call readi(mcmcard,'NSEED',nseed,20)
3417 call readi(mcmcard,'NTOTAL',ntotal,10000)
3418 call reada(mcmcard,'CUT1',cut1,2.0d0)
3419 call reada(mcmcard,'CUT2',cut2,5.0d0)
3420 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3421 call readi(mcmcard,'ICMAX',icmax,3)
3422 call readi(mcmcard,'IRESTART',irestart,0)
3423 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3426 call reada(mcmcard,'DELE',dele,20.0d0)
3427 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3428 call readi(mcmcard,'IREF',iref,0)
3429 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3430 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3431 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3432 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3433 write (iout,*) "NCONF_IN",nconf_in
3435 end subroutine csaread
3436 !-----------------------------------------------------------------------------
3440 use control_data, only: MaxMoveType
3443 ! implicit real*8 (a-h,o-z)
3444 ! include 'DIMENSIONS'
3445 ! include 'COMMON.MCM'
3446 ! include 'COMMON.MCE'
3447 ! include 'COMMON.IOUNITS'
3448 ! character(len=80) :: ucase
3449 character(len=320) :: mcmcard
3452 ! real(kind=8) :: var,ene
3454 call card_concat(mcmcard,.true.)
3455 call readi(mcmcard,'MAXACC',maxacc,100)
3456 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3457 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3458 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3459 call readi(mcmcard,'MAXREPM',maxrepm,200)
3460 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3461 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3462 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3463 call reada(mcmcard,'E_UP',e_up,5.0D0)
3464 call reada(mcmcard,'DELTE',delte,0.1D0)
3465 call readi(mcmcard,'NSWEEP',nsweep,5)
3466 call readi(mcmcard,'NSTEPH',nsteph,0)
3467 call readi(mcmcard,'NSTEPC',nstepc,0)
3468 call reada(mcmcard,'TMIN',tmin,298.0D0)
3469 call reada(mcmcard,'TMAX',tmax,298.0D0)
3470 call readi(mcmcard,'NWINDOW',nwindow,0)
3471 call readi(mcmcard,'PRINT_MC',print_mc,0)
3472 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3473 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3474 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3475 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3476 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3477 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3478 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3479 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3480 if (nwindow.gt.0) then
3481 allocate(winstart(nwindow)) !!el (maxres)
3482 allocate(winend(nwindow)) !!el
3483 allocate(winlen(nwindow)) !!el
3484 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3486 winlen(i)=winend(i)-winstart(i)+1
3489 if (tmax.lt.tmin) tmax=tmin
3490 if (tmax.eq.tmin) then
3494 if (nstepc.gt.0 .and. nsteph.gt.0) then
3495 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3496 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3498 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3499 ! Probabilities of different move types
3500 sumpro_type(0)=0.0D0
3501 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3502 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3503 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3504 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3505 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3506 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3507 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3509 print *,'i',i,' sumprotype',sumpro_type(i)
3510 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3511 print *,'i',i,' sumprotype',sumpro_type(i)
3514 end subroutine mcmread
3515 !-----------------------------------------------------------------------------
3516 subroutine read_minim
3519 ! implicit real*8 (a-h,o-z)
3520 ! include 'DIMENSIONS'
3521 ! include 'COMMON.MINIM'
3522 ! include 'COMMON.IOUNITS'
3523 ! character(len=80) :: ucase
3524 character(len=320) :: minimcard
3526 ! integer :: ntf,ik,iw_pdb
3527 ! real(kind=8) :: var,ene
3529 call card_concat(minimcard,.true.)
3530 call readi(minimcard,'MAXMIN',maxmin,2000)
3531 call readi(minimcard,'MAXFUN',maxfun,5000)
3532 call readi(minimcard,'MINMIN',minmin,maxmin)
3533 call readi(minimcard,'MINFUN',minfun,maxmin)
3534 call reada(minimcard,'TOLF',tolf,1.0D-2)
3535 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3536 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3537 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3538 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3539 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3540 'Options in energy minimization:'
3541 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3542 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3543 'MinMin:',MinMin,' MinFun:',MinFun,&
3544 ' TolF:',TolF,' RTolF:',RTolF
3546 end subroutine read_minim
3547 !-----------------------------------------------------------------------------
3548 subroutine openunits
3550 use MD_data, only: usampl
3553 use control_data, only:out1file
3554 use control, only: getenv_loc
3555 ! implicit real*8 (a-h,o-z)
3556 ! include 'DIMENSIONS'
3559 character(len=16) :: form,nodename
3560 integer :: nodelen,ierror,npos
3562 ! include 'COMMON.SETUP'
3563 ! include 'COMMON.IOUNITS'
3564 ! include 'COMMON.MD'
3565 ! include 'COMMON.CONTROL'
3566 integer :: lenpre,lenpot,lentmp !,ilen
3568 character(len=3) :: out1file_text !,ucase
3569 character(len=3) :: ll
3572 ! integer :: ntf,ik,iw_pdb
3573 ! real(kind=8) :: var,ene
3575 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3576 call getenv_loc("PREFIX",prefix)
3578 call getenv_loc("POT",pot)
3579 call getenv_loc("DIRTMP",tmpdir)
3580 call getenv_loc("CURDIR",curdir)
3581 call getenv_loc("OUT1FILE",out1file_text)
3582 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3583 out1file_text=ucase(out1file_text)
3584 if (out1file_text(1:1).eq."Y") then
3587 out1file=fg_rank.gt.0
3592 if (lentmp.gt.0) then
3593 write (*,'(80(1h!))')
3594 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3595 write (*,'(80(1h!))')
3596 write (*,*)"All output files will be on node /tmp directory."
3598 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3599 if (me.eq.king) then
3600 write (*,*) "The master node is ",nodename
3601 else if (fg_rank.eq.0) then
3602 write (*,*) "I am the CG slave node ",nodename
3604 write (*,*) "I am the FG slave node ",nodename
3607 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3608 lenpre = lentmp+lenpre+1
3610 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3611 ! Get the names and open the input files
3612 #if defined(WINIFL) || defined(WINPGI)
3613 open(1,file=pref_orig(:ilen(pref_orig))// &
3614 '.inp',status='old',readonly,shared)
3615 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3616 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3617 ! Get parameter filenames and open the parameter files.
3618 call getenv_loc('BONDPAR',bondname)
3619 open (ibond,file=bondname,status='old',readonly,shared)
3620 call getenv_loc('THETPAR',thetname)
3621 open (ithep,file=thetname,status='old',readonly,shared)
3622 call getenv_loc('ROTPAR',rotname)
3623 open (irotam,file=rotname,status='old',readonly,shared)
3624 call getenv_loc('TORPAR',torname)
3625 open (itorp,file=torname,status='old',readonly,shared)
3626 call getenv_loc('TORDPAR',tordname)
3627 open (itordp,file=tordname,status='old',readonly,shared)
3628 call getenv_loc('FOURIER',fouriername)
3629 open (ifourier,file=fouriername,status='old',readonly,shared)
3630 call getenv_loc('ELEPAR',elename)
3631 open (ielep,file=elename,status='old',readonly,shared)
3632 call getenv_loc('SIDEPAR',sidename)
3633 open (isidep,file=sidename,status='old',readonly,shared)
3634 #elif (defined CRAY) || (defined AIX)
3635 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3637 ! print *,"Processor",myrank," opened file 1"
3638 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3639 ! print *,"Processor",myrank," opened file 9"
3640 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3641 ! Get parameter filenames and open the parameter files.
3642 call getenv_loc('BONDPAR',bondname)
3643 open (ibond,file=bondname,status='old',action='read')
3644 ! print *,"Processor",myrank," opened file IBOND"
3645 call getenv_loc('THETPAR',thetname)
3646 open (ithep,file=thetname,status='old',action='read')
3647 ! print *,"Processor",myrank," opened file ITHEP"
3648 call getenv_loc('ROTPAR',rotname)
3649 open (irotam,file=rotname,status='old',action='read')
3650 ! print *,"Processor",myrank," opened file IROTAM"
3651 call getenv_loc('TORPAR',torname)
3652 open (itorp,file=torname,status='old',action='read')
3653 ! print *,"Processor",myrank," opened file ITORP"
3654 call getenv_loc('TORDPAR',tordname)
3655 open (itordp,file=tordname,status='old',action='read')
3656 ! print *,"Processor",myrank," opened file ITORDP"
3657 call getenv_loc('SCCORPAR',sccorname)
3658 open (isccor,file=sccorname,status='old',action='read')
3659 ! print *,"Processor",myrank," opened file ISCCOR"
3660 call getenv_loc('FOURIER',fouriername)
3661 open (ifourier,file=fouriername,status='old',action='read')
3662 ! print *,"Processor",myrank," opened file IFOURIER"
3663 call getenv_loc('ELEPAR',elename)
3664 open (ielep,file=elename,status='old',action='read')
3665 ! print *,"Processor",myrank," opened file IELEP"
3666 call getenv_loc('SIDEPAR',sidename)
3667 open (isidep,file=sidename,status='old',action='read')
3668 ! print *,"Processor",myrank," opened file ISIDEP"
3669 ! print *,"Processor",myrank," opened parameter files"
3671 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3672 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3673 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3674 ! Get parameter filenames and open the parameter files.
3675 call getenv_loc('BONDPAR',bondname)
3676 open (ibond,file=bondname,status='old')
3677 call getenv_loc('THETPAR',thetname)
3678 open (ithep,file=thetname,status='old')
3679 call getenv_loc('ROTPAR',rotname)
3680 open (irotam,file=rotname,status='old')
3681 call getenv_loc('TORPAR',torname)
3682 open (itorp,file=torname,status='old')
3683 call getenv_loc('TORDPAR',tordname)
3684 open (itordp,file=tordname,status='old')
3685 call getenv_loc('SCCORPAR',sccorname)
3686 open (isccor,file=sccorname,status='old')
3687 call getenv_loc('FOURIER',fouriername)
3688 open (ifourier,file=fouriername,status='old')
3689 call getenv_loc('ELEPAR',elename)
3690 open (ielep,file=elename,status='old')
3691 call getenv_loc('SIDEPAR',sidename)
3692 open (isidep,file=sidename,status='old')
3694 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3696 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3697 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3698 ! Get parameter filenames and open the parameter files.
3699 call getenv_loc('BONDPAR',bondname)
3700 open (ibond,file=bondname,status='old',action='read')
3701 call getenv_loc('THETPAR',thetname)
3702 open (ithep,file=thetname,status='old',action='read')
3703 call getenv_loc('ROTPAR',rotname)
3704 open (irotam,file=rotname,status='old',action='read')
3705 call getenv_loc('TORPAR',torname)
3706 open (itorp,file=torname,status='old',action='read')
3707 call getenv_loc('TORDPAR',tordname)
3708 open (itordp,file=tordname,status='old',action='read')
3709 call getenv_loc('SCCORPAR',sccorname)
3710 open (isccor,file=sccorname,status='old',action='read')
3712 call getenv_loc('THETPARPDB',thetname_pdb)
3713 print *,"thetname_pdb ",thetname_pdb
3714 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3715 print *,ithep_pdb," opened"
3717 call getenv_loc('FOURIER',fouriername)
3718 open (ifourier,file=fouriername,status='old',readonly)
3719 call getenv_loc('ELEPAR',elename)
3720 open (ielep,file=elename,status='old',readonly)
3721 call getenv_loc('SIDEPAR',sidename)
3722 open (isidep,file=sidename,status='old',readonly)
3724 call getenv_loc('ROTPARPDB',rotname_pdb)
3725 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3730 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3731 ! Use -DOLDSCP to use hard-coded constants instead.
3733 call getenv_loc('SCPPAR',scpname)
3734 #if defined(WINIFL) || defined(WINPGI)
3735 open (iscpp,file=scpname,status='old',readonly,shared)
3736 #elif (defined CRAY) || (defined AIX)
3737 open (iscpp,file=scpname,status='old',action='read')
3739 open (iscpp,file=scpname,status='old')
3741 open (iscpp,file=scpname,status='old',action='read')
3744 call getenv_loc('PATTERN',patname)
3745 #if defined(WINIFL) || defined(WINPGI)
3746 open (icbase,file=patname,status='old',readonly,shared)
3747 #elif (defined CRAY) || (defined AIX)
3748 open (icbase,file=patname,status='old',action='read')
3750 open (icbase,file=patname,status='old')
3752 open (icbase,file=patname,status='old',action='read')
3755 ! Open output file only for CG processes
3756 ! print *,"Processor",myrank," fg_rank",fg_rank
3757 if (fg_rank.eq.0) then
3759 if (nodes.eq.1) then
3762 npos = dlog10(dfloat(nodes-1))+1
3764 if (npos.lt.3) npos=3
3765 write (liczba,'(i1)') npos
3766 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3768 write (liczba,form) me
3769 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3770 liczba(:ilen(liczba))
3771 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3773 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3775 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3776 liczba(:ilen(liczba))//'.mol2'
3777 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3778 liczba(:ilen(liczba))//'.stat'
3780 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3781 //liczba(:ilen(liczba))//'.stat')
3782 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3785 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3786 liczba(:ilen(liczba))//'.const'
3791 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3792 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3793 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3794 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3795 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3797 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3799 rest2name=prefix(:ilen(prefix))//'.rst'
3801 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3804 #if defined(AIX) || defined(PGI)
3805 if (me.eq.king .or. .not. out1file) &
3806 open(iout,file=outname,status='unknown')
3808 if (fg_rank.gt.0) then
3809 write (liczba,'(i3.3)') myrank/nfgtasks
3810 write (ll,'(bz,i3.3)') fg_rank
3811 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3816 open(igeom,file=intname,status='unknown',position='append')
3817 open(ipdb,file=pdbname,status='unknown')
3818 open(imol2,file=mol2name,status='unknown')
3819 open(istat,file=statname,status='unknown',position='append')
3821 !1out open(iout,file=outname,status='unknown')
3824 if (me.eq.king .or. .not.out1file) &
3825 open(iout,file=outname,status='unknown')
3827 if (fg_rank.gt.0) then
3828 write (liczba,'(i3.3)') myrank/nfgtasks
3829 write (ll,'(bz,i3.3)') fg_rank
3830 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3835 open(igeom,file=intname,status='unknown',access='append')
3836 open(ipdb,file=pdbname,status='unknown')
3837 open(imol2,file=mol2name,status='unknown')
3838 open(istat,file=statname,status='unknown',access='append')
3840 !1out open(iout,file=outname,status='unknown')
3843 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3844 csa_seed=prefix(:lenpre)//'.CSA.seed'
3845 csa_history=prefix(:lenpre)//'.CSA.history'
3846 csa_bank=prefix(:lenpre)//'.CSA.bank'
3847 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3848 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3849 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3850 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3851 csa_int=prefix(:lenpre)//'.int'
3852 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3853 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3854 csa_in=prefix(:lenpre)//'.CSA.in'
3855 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3858 write (iout,'(80(1h-))')
3859 write (iout,'(30x,a)') "FILE ASSIGNMENT"
3860 write (iout,'(80(1h-))')
3861 write (iout,*) "Input file : ",&
3862 pref_orig(:ilen(pref_orig))//'.inp'
3863 write (iout,*) "Output file : ",&
3864 outname(:ilen(outname))
3866 write (iout,*) "Sidechain potential file : ",&
3867 sidename(:ilen(sidename))
3869 write (iout,*) "SCp potential file : ",&
3870 scpname(:ilen(scpname))
3872 write (iout,*) "Electrostatic potential file : ",&
3873 elename(:ilen(elename))
3874 write (iout,*) "Cumulant coefficient file : ",&
3875 fouriername(:ilen(fouriername))
3876 write (iout,*) "Torsional parameter file : ",&
3877 torname(:ilen(torname))
3878 write (iout,*) "Double torsional parameter file : ",&
3879 tordname(:ilen(tordname))
3880 write (iout,*) "SCCOR parameter file : ",&
3881 sccorname(:ilen(sccorname))
3882 write (iout,*) "Bond & inertia constant file : ",&
3883 bondname(:ilen(bondname))
3884 write (iout,*) "Bending parameter file : ",&
3885 thetname(:ilen(thetname))
3886 write (iout,*) "Rotamer parameter file : ",&
3887 rotname(:ilen(rotname))
3890 write (iout,*) "Thetpdb parameter file : ",&
3891 thetname_pdb(:ilen(thetname_pdb))
3894 write (iout,*) "Threading database : ",&
3895 patname(:ilen(patname))
3897 write (iout,*)" DIRTMP : ",&
3899 write (iout,'(80(1h-))')
3902 end subroutine openunits
3903 !-----------------------------------------------------------------------------
3906 use geometry_data, only: nres,dc
3908 ! implicit real*8 (a-h,o-z)
3909 ! include 'DIMENSIONS'
3910 ! include 'COMMON.CHAIN'
3911 ! include 'COMMON.IOUNITS'
3912 ! include 'COMMON.MD'
3915 ! real(kind=8) :: var,ene
3917 open(irest2,file=rest2name,status='unknown')
3918 read(irest2,*) totT,EK,potE,totE,t_bath
3920 ! AL 4/17/17: Now reading d_t(0,:) too
3922 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
3925 ! AL 4/17/17: Now reading d_c(0,:) too
3927 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
3930 read (irest2,*) iset
3934 end subroutine readrst
3935 !-----------------------------------------------------------------------------
3936 subroutine read_fragments
3940 use control_data, only:out1file
3943 ! implicit real*8 (a-h,o-z)
3944 ! include 'DIMENSIONS'
3948 ! include 'COMMON.SETUP'
3949 ! include 'COMMON.CHAIN'
3950 ! include 'COMMON.IOUNITS'
3951 ! include 'COMMON.MD'
3952 ! include 'COMMON.CONTROL'
3955 ! real(kind=8) :: var,ene
3957 read(inp,*) nset,nfrag,npair,nfrag_back
3959 !el from module energy
3960 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
3961 if(.not.allocated(wfrag_back)) then
3962 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3963 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3965 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
3966 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
3968 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
3969 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
3972 if(me.eq.king.or..not.out1file) &
3973 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
3974 " nfrag_back",nfrag_back
3976 read(inp,*) mset(iset)
3978 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
3980 if(me.eq.king.or..not.out1file) &
3981 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
3982 ifrag(2,i,iset), qinfrag(i,iset)
3985 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
3987 if(me.eq.king.or..not.out1file) &
3988 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
3989 ipair(2,i,iset), qinpair(i,iset)
3992 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3993 wfrag_back(3,i,iset),&
3994 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
3995 if(me.eq.king.or..not.out1file) &
3996 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3997 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4001 end subroutine read_fragments
4002 !-----------------------------------------------------------------------------
4004 !-----------------------------------------------------------------------------
4008 ! implicit real*8 (a-h,o-z)
4009 ! include 'DIMENSIONS'
4010 ! include 'COMMON.CSA'
4011 ! include 'COMMON.BANK'
4012 ! include 'COMMON.IOUNITS'
4014 ! integer :: ntf,ik,iw_pdb
4015 ! real(kind=8) :: var,ene
4017 open(icsa_in,file=csa_in,status="old",err=100)
4018 read(icsa_in,*) nconf
4019 read(icsa_in,*) jstart,jend
4020 read(icsa_in,*) nstmax
4021 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4022 read(icsa_in,*) nran0,nran1,irr
4023 read(icsa_in,*) nseed
4024 read(icsa_in,*) ntotal,cut1,cut2
4025 read(icsa_in,*) estop
4026 read(icsa_in,*) icmax,irestart
4027 read(icsa_in,*) ntbankm,dele,difcut
4028 read(icsa_in,*) iref,rmscut,pnccut
4029 read(icsa_in,*) ndiff
4036 end subroutine csa_read
4037 !-----------------------------------------------------------------------------
4038 subroutine initial_write
4041 ! implicit real*8 (a-h,o-z)
4042 ! include 'DIMENSIONS'
4043 ! include 'COMMON.CSA'
4044 ! include 'COMMON.BANK'
4045 ! include 'COMMON.IOUNITS'
4047 ! integer :: ntf,ik,iw_pdb
4048 ! real(kind=8) :: var,ene
4050 open(icsa_seed,file=csa_seed,status="unknown")
4051 write(icsa_seed,*) "seed"
4053 #if defined(AIX) || defined(PGI)
4054 open(icsa_history,file=csa_history,status="unknown",&
4057 open(icsa_history,file=csa_history,status="unknown",&
4060 write(icsa_history,*) nconf
4061 write(icsa_history,*) jstart,jend
4062 write(icsa_history,*) nstmax
4063 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4064 write(icsa_history,*) nran0,nran1,irr
4065 write(icsa_history,*) nseed
4066 write(icsa_history,*) ntotal,cut1,cut2
4067 write(icsa_history,*) estop
4068 write(icsa_history,*) icmax,irestart
4069 write(icsa_history,*) ntbankm,dele,difcut
4070 write(icsa_history,*) iref,rmscut,pnccut
4071 write(icsa_history,*) ndiff
4073 write(icsa_history,*)
4076 open(icsa_bank1,file=csa_bank1,status="unknown")
4077 write(icsa_bank1,*) 0
4081 end subroutine initial_write
4082 !-----------------------------------------------------------------------------
4083 subroutine restart_write
4086 ! implicit real*8 (a-h,o-z)
4087 ! include 'DIMENSIONS'
4088 ! include 'COMMON.IOUNITS'
4089 ! include 'COMMON.CSA'
4090 ! include 'COMMON.BANK'
4092 ! integer :: ntf,ik,iw_pdb
4093 ! real(kind=8) :: var,ene
4095 #if defined(AIX) || defined(PGI)
4096 open(icsa_history,file=csa_history,position="append")
4098 open(icsa_history,file=csa_history,access="append")
4100 write(icsa_history,*)
4101 write(icsa_history,*) "This is restart"
4102 write(icsa_history,*)
4103 write(icsa_history,*) nconf
4104 write(icsa_history,*) jstart,jend
4105 write(icsa_history,*) nstmax
4106 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4107 write(icsa_history,*) nran0,nran1,irr
4108 write(icsa_history,*) nseed
4109 write(icsa_history,*) ntotal,cut1,cut2
4110 write(icsa_history,*) estop
4111 write(icsa_history,*) icmax,irestart
4112 write(icsa_history,*) ntbankm,dele,difcut
4113 write(icsa_history,*) iref,rmscut,pnccut
4114 write(icsa_history,*) ndiff
4115 write(icsa_history,*)
4116 write(icsa_history,*) "irestart is: ", irestart
4118 write(icsa_history,*)
4122 end subroutine restart_write
4123 !-----------------------------------------------------------------------------
4125 !-----------------------------------------------------------------------------
4126 subroutine write_pdb(npdb,titelloc,ee)
4128 ! implicit real*8 (a-h,o-z)
4129 ! include 'DIMENSIONS'
4130 ! include 'COMMON.IOUNITS'
4131 character(len=50) :: titelloc1
4132 character*(*) :: titelloc
4133 character(len=3) :: zahl
4134 character(len=5) :: liczba5
4136 integer :: npdb !,ilen
4140 ! real(kind=8) :: var,ene
4144 if (npdb.lt.1000) then
4145 call numstr(npdb,zahl)
4146 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4148 if (npdb.lt.10000) then
4149 write(liczba5,'(i1,i4)') 0,npdb
4151 write(liczba5,'(i5)') npdb
4153 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4155 call pdbout(ee,titelloc1,ipdb)
4158 end subroutine write_pdb
4159 !-----------------------------------------------------------------------------
4161 !-----------------------------------------------------------------------------
4162 subroutine write_thread_summary
4163 ! Thread the sequence through a database of known structures
4164 use control_data, only: refstr
4166 use energy_data, only: n_ene_comp
4168 ! implicit real*8 (a-h,o-z)
4169 ! include 'DIMENSIONS'
4171 use MPI_data !include 'COMMON.INFO'
4174 ! include 'COMMON.CONTROL'
4175 ! include 'COMMON.CHAIN'
4176 ! include 'COMMON.DBASE'
4177 ! include 'COMMON.INTERACT'
4178 ! include 'COMMON.VAR'
4179 ! include 'COMMON.THREAD'
4180 ! include 'COMMON.FFIELD'
4181 ! include 'COMMON.SBRIDGE'
4182 ! include 'COMMON.HEADER'
4183 ! include 'COMMON.NAMES'
4184 ! include 'COMMON.IOUNITS'
4185 ! include 'COMMON.TIME1'
4187 integer,dimension(maxthread) :: ip
4188 real(kind=8),dimension(0:n_ene) :: energia
4190 integer :: i,j,ii,jj,ipj,ik,kk,ist
4191 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4193 write (iout,'(30x,a/)') &
4194 ' *********** Summary threading statistics ************'
4195 write (iout,'(a)') 'Initial energies:'
4196 write (iout,'(a4,2x,a12,14a14,3a8)') &
4197 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4198 'RMSnat','NatCONT','NNCONT','RMS'
4199 ! Energy sort patterns
4204 enet=ener(n_ene-1,ip(i))
4207 if (ener(n_ene-1,ip(j)).lt.enet) then
4209 enet=ener(n_ene-1,ip(j))
4221 ist=nres_base(2,ii)+ipatt(2,i)
4223 energia(i)=ener0(kk,i)
4225 etot=ener0(n_ene_comp+1,i)
4226 rmsnat=ener0(n_ene_comp+2,i)
4227 rms=ener0(n_ene_comp+3,i)
4228 frac=ener0(n_ene_comp+4,i)
4229 frac_nn=ener0(n_ene_comp+5,i)
4232 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4233 i,str_nam(ii),ist+1,&
4234 (energia(print_order(kk)),kk=1,nprint_ene),&
4235 etot,rmsnat,frac,frac_nn,rms
4237 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4238 i,str_nam(ii),ist+1,&
4239 (energia(print_order(kk)),kk=1,nprint_ene),etot
4242 write (iout,'(//a)') 'Final energies:'
4243 write (iout,'(a4,2x,a12,17a14,3a8)') &
4244 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4245 'RMSnat','NatCONT','NNCONT','RMS'
4249 ist=nres_base(2,ii)+ipatt(2,i)
4251 energia(kk)=ener(kk,ik)
4253 etot=ener(n_ene_comp+1,i)
4254 rmsnat=ener(n_ene_comp+2,i)
4255 rms=ener(n_ene_comp+3,i)
4256 frac=ener(n_ene_comp+4,i)
4257 frac_nn=ener(n_ene_comp+5,i)
4258 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4259 i,str_nam(ii),ist+1,&
4260 (energia(print_order(kk)),kk=1,nprint_ene),&
4261 etot,rmsnat,frac,frac_nn,rms
4263 write (iout,'(/a/)') 'IEXAM array:'
4264 write (iout,'(i5)') nexcl
4266 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4268 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4269 'Max. time for threading step ',max_time_for_thread,&
4270 'Average time for threading step: ',ave_time_for_thread
4272 end subroutine write_thread_summary
4273 !-----------------------------------------------------------------------------
4274 subroutine write_stat_thread(ithread,ipattern,ist)
4276 use energy_data, only: n_ene_comp
4278 ! implicit real*8 (a-h,o-z)
4279 ! include "DIMENSIONS"
4280 ! include "COMMON.CONTROL"
4281 ! include "COMMON.IOUNITS"
4282 ! include "COMMON.THREAD"
4283 ! include "COMMON.FFIELD"
4284 ! include "COMMON.DBASE"
4285 ! include "COMMON.NAMES"
4286 real(kind=8),dimension(0:n_ene) :: energia
4288 integer :: ithread,ipattern,ist,i
4289 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4291 #if defined(AIX) || defined(PGI)
4292 open(istat,file=statname,position='append')
4294 open(istat,file=statname,access='append')
4297 energia(i)=ener(i,ithread)
4299 etot=ener(n_ene_comp+1,ithread)
4300 rmsnat=ener(n_ene_comp+2,ithread)
4301 rms=ener(n_ene_comp+3,ithread)
4302 frac=ener(n_ene_comp+4,ithread)
4303 frac_nn=ener(n_ene_comp+5,ithread)
4304 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4305 ithread,str_nam(ipattern),ist+1,&
4306 (energia(print_order(i)),i=1,nprint_ene),&
4307 etot,rmsnat,frac,frac_nn,rms
4310 end subroutine write_stat_thread
4311 !-----------------------------------------------------------------------------
4313 !-----------------------------------------------------------------------------
4314 end module io_config